text:
# Lepton-flavor-violating decays 𝐿→𝑙𝛾𝛾 as a new probe of supersymmetry with broken R-parity

thanks: Work supported in part by the EC 5th framework, contract number HPMF-CT-2002-01663

## Abstract

Lepton-flavor-violating decays of the type $`Ll\gamma \gamma `$ ($`\mu e\gamma \gamma `$, $`\tau e\gamma \gamma `$, and $`\tau \mu \gamma \gamma `$) are proposed as new probes of R-parity-violating supersymmetry. Non-penguin diagrams with a sneutrino that decays into two photons via a triangle graph might trigger such decays even in the absence of the corresponding radiative decays into one photon only, e.g. $`\mu e\gamma `$. Thus, processes of the type $`Ll\gamma \gamma `$ may provide an independent probe of new flavor physics.

preprint: CERN-TH/2003-034

Supersymmetry (SUSY) is one of the most promising candidates for new physics, curing some of the shortcomings of the Standard Model (SM), for instance, the emergence of quadratic divergences in the Higgs sector SUSY . Despite its appealing theoretical features, however, there is so far no experimental evidence for SUSY up to the electroweak scale. Apart from direct evidence, i.e. direct production and decays of SUSY particles, indirect probes could be employed to search for SUSY, i.e. processes in which the SUSY particles emerge virtually either in loops or as tree-level mediators. In this respect Flavor-Changing Neutral Current (FCNC) or Lepton-Flavor-Violating (LFV) processes are a natural searching ground for SUSY, since the SM background to such processes is either loop-suppressed (FCNC) or does not exist (LFV).

In this paper the effects of new physics on a class of LFV processes, $`Ll\gamma \gamma `$, are investigated, in which a lepton decays radiatively into a lighter one plus two photons. As an example, R-Parity-Violating (RPV) SUSY is discussed.

In the conventional RPV framework, the SUSY superpotential is supplemented with new RPV interaction terms that are either baryon- or lepton-number-violating Drei ; RPV . For the purpose of this paper, only the following lepton-number-violating operators are relevant:

$`𝒲_{RPV}{\displaystyle \frac{1}{2}}\lambda _{ijk}ϵ_{ab}\widehat{L}_i^a\widehat{L}_j^b\widehat{E}_k^c+\lambda _{ijk}^{}ϵ_{ab}\widehat{L}_i^a\widehat{Q}_j^b\widehat{D}_k^c,`$ (1)
where $`\widehat{Q}`$ and $`\widehat{L}`$ are SU(2) doublet quark and lepton supermultiplets, respectively, and $`\widehat{D}^c`$ and $`\widehat{E}^c`$ denote the SU(2) singlet down-type quark and lepton supermultiplet; $`i,j,k=1,2,3`$ label the generations and $`a,b`$ are the SU(2) indices that are responsible for the antisymmetric property $`\lambda _{ijk}=\lambda _{jik}`$.

Within the framework of RPV SUSY, $`\mu e\gamma `$ (and $`\tau e\gamma `$, $`\tau \mu \gamma `$) can be mediated by penguin-like diagrams of the type depicted in Fig. 1.

Depending on the RPV couplings involved, the particles running in the loop are either quark–squark, or lepton–slepton. The process $`Ll\gamma `$, where $`L`$ and $`l`$ are leptons of generations $`i`$ and $`j`$, respectively, is then proportional to either $`\lambda _{imn}\lambda _{jmn}`$ or $`\lambda _{mni}\lambda _{mnj}`$ or $`\lambda _{imn}^{}\lambda _{jmn}^{}`$ Chai .

So far, very little attention has been devoted to $`Ll\gamma \gamma `$, under the assumption that it is always suppressed with respect to $`Ll\gamma `$. An exception to this rule is discussed in leptoq , where a leptoquark model in which $`Ll\gamma \gamma `$ is much larger than $`Ll\gamma `$ is presented. Note that their model yields $`BR(\mu e\gamma \gamma )<10^{18}`$; see also Bowman . In general, an additional photon can easily be accommodated - one just has to add a photon to the diagram in Fig. 1 and obtain the $`\lambda `$-reducible diagrams shown in Fig. 2a (they are reducible in the sense that, after removing a photon line from $`Ll\gamma \gamma `$, a legitimate RPV diagram for $`Ll\gamma `$ is obtained). Obviously processes mediated through such diagrams do not provide any further information on the corresponding RPV couplings, since they differ only by kinematic factors and by an extra electromagnetic coupling.

However, it is quite possible that the decay $`Ll\gamma `$ does not exist, or its rate is too small to be observed, whereas the decay $`Ll\gamma \gamma `$ still has an appreciable rate. This could be the case if, for instance, all products of the type $`\lambda \lambda `$ and $`\lambda ^{}\lambda ^{}`$ are much smaller than products of the type $`\lambda \lambda ^{}`$, in which case the $`\lambda `$-irreducible topology of Fig. 2b may give rise to a much larger rate than the $`\lambda `$-reducible diagrams of Fig. 2a.

For example, if $`\lambda _{122}`$ and $`\lambda _{233}^{}`$ are the only non-vanishing lambdas, then $`\mathrm{\Gamma }(\mu e\gamma )=0`$ while $`\mathrm{\Gamma }(\mu e\gamma \gamma )0`$. Accordingly, such RPV coupling product combinations can be detected or constrained by the measured upper limit on the branching ratio of $`\mu e\gamma \gamma `$, which is at present $`7.2\times 10^{11}`$ PDG ; experiment2 . No experimental limits exist so far for the analogous decays of the $`\tau `$-lepton.

There are proposals to measure (or constrain) the branching ratios of such SM-forbidden decays of the $`\mu `$ down to branching ratios of the order of $`10^{14}`$ experiment4 , $`10^{15}`$ experiment5  or even down to $`10^{18}`$ experiment6 . This highlights the phenomenological impact of LFV muon decays.

A general basis of effective operators describing the decays $`Ll\gamma \gamma `$ with the topology of Fig. 2, i.e. $`LlS_kl\gamma \gamma `$, will be constructed first. Here $`S_k`$ is a scalar particle of type $`k`$ and $`i,j`$ are generation indices of $`L,l`$, respectively. The momentum notation $`L(p)l(p^{})\gamma (k)\gamma (k^{})`$ will be used.

The amplitude can be defined as

$$_{Ll\gamma \gamma }=\underset{n=1}{\overset{2}{}}\underset{𝒫=L,R}{}C_{n(ijkm)}^𝒫𝒪_n^𝒫,$$
(2)
where $`m`$ is the family index of the fermion that runs in the loop, and the index $`k`$ corresponds to the intermediate scalar. Also, $`C_{n(ijkm)}^L,C_{n(ijkm)}^R`$ are Wilson-like coefficients with $`L(R)=(1(+)\gamma _5)/2`$, and

$`𝒪_1^L`$ $`=`$ $`\left[\overline{u}_l(p^{})Lu_L(p)\right]\left(kk^{}g^{\lambda \nu }k^\lambda k^\nu \right)ϵ_\lambda ^{}ϵ_\nu ^{}`$
$`𝒪_2^L`$ $`=`$ $`\left[\overline{u}_l(p^{})Lu_L(p)\right]i\epsilon ^{\lambda \nu \rho \sigma }k_\rho k_\sigma ^{}ϵ_\lambda ^{}ϵ_\nu ^{}`$ (3)
are effective operators. If the incoming lepton is right-handed, then $`LR`$ in (Lepton-flavor-violating decays $`Ll\gamma \gamma `$ as a new probe of supersymmetry with broken R-parity).

In the framework of RPV SUSY with a sneutrino ($`\stackrel{~}{\nu }_k`$) exchange, the Wilson coefficients are given by:

$`C_{1;2(ijkm)}^L;C_{1;2(ijkm)}^R(\mathrm{}_m\mathrm{in}\mathrm{the}\mathrm{loop})`$ $`=`$ $`0;{\displaystyle \frac{\alpha }{4\pi }}{\displaystyle \frac{i\lambda _{kji}^{}\lambda _{kmm}}{m_\mathrm{}_mM_{\stackrel{~}{\nu }_k}^2}}[f_{1/2}(x);g_{1/2}(x)],`$
$`C_{1;2(ijkm)}^L;C_{1;2(ijkm)}^R(d_m\mathrm{in}\mathrm{the}\mathrm{loop})`$ $`=`$ $`0;N_cQ_d^2{\displaystyle \frac{\alpha }{4\pi }}{\displaystyle \frac{i\lambda _{kji}^{}\lambda _{kmm}^{}}{m_{d_m}M_{\stackrel{~}{\nu }_k}^2}}[f_{1/2}(x);g_{1/2}(x)],`$ (4)
where $`M_{\stackrel{~}{\nu }_k}`$ is the mass of the $`k`$th generation sneutrino and $`m_f`$, $`f=\mathrm{}\mathrm{or}d`$, is the mass of the fermion of generation $`m`$ that runs in the loop. Furthermore $`N_c=3`$, $`Q_d=1/3`$ and the abbreviation $`x=2m_f^2/kk^{}`$ is being used. For diagrams with a charge-conjugate sneutrino ($`\stackrel{~}{\nu }^{}`$) exchange, in (4) $`C_{1;2(ijkm)}^LC_{1;2(jikm)}^R`$ (i.e. with $`\lambda _{kji}^{}\lambda _{kij}^{}`$) and $`M_{\stackrel{~}{\nu }_k}M_{\stackrel{~}{\nu }_k^{}}`$. The functions $`f_{1/2}`$ and $`g_{1/2}`$ (originating from the loop integral of the triangle) are given by

$`f_{1/2}(x)`$ $`=`$ $`2x\left[1+(1x)\mathrm{arcsin}^2\left({\displaystyle \frac{1}{\sqrt{x}}}\right)\right],`$
$`g_{1/2}(x)`$ $`=`$ $`2x\mathrm{arcsin}^2\left({\displaystyle \frac{1}{\sqrt{x}}}\right).`$ (5)
The emergence of these functions reflects the different structure of the integral related to the parity of the couplings of the fermion in the loop. The subscripts refer to the spin of the particle running in the triangle. Note that when $`x<1`$, the $`\mathrm{arcsin}\left(1/\sqrt{x}\right)`$ develops an imaginary part, which corresponds to the loop particle residing on its mass shell.

In what follows, results will be given for both constituent and current light quark masses: $`m_d=300`$ MeV, $`m_s=450`$ MeV, $`m_b=4.5`$ GeV, and $`m_d=10`$ MeV, $`m_s=120`$ MeV, $`m_b=4.2`$ GeV, respectively. We believe that using constituent quark masses makes more sense, since for our case here, it would allow the decay $`\mu ed\overline{d}`$, for current $`d`$ quarks Eeg .

For the case of the muon decay, $`\mu e\gamma \gamma `$, and for either down quarks with constituent masses or for the $`\tau `$-lepton running in the loop, the relation $`m_f^2>m_\mu ^2>kk^{}`$ holds. Then, the limit $`m_f^2kk^{}`$ ($`m_f=m_d,m_s,m_b`$ or $`m_\tau `$) provides a good approximation. Therefore, the loop integral functions can be replaced with their $`x\mathrm{}`$ limits, $`f_{1/2}(x\mathrm{})=\frac{4}{3}`$ and $`g_{1/2}(x\mathrm{})=2`$. Using this approximation, decay widths are given by

$`\mathrm{\Gamma }(\mu e\gamma \gamma )_{m_f^2kk^{}}A_{RPV}^{fS}\left({\displaystyle \frac{\alpha }{\pi }}\right)^2{\displaystyle \frac{1}{256\pi ^3}}{\displaystyle \frac{13m_\mu ^7}{4320m_f^2M_S^4}},`$ (6)
where $`f=d_i`$ ($`i=1,2,3`$ for $`d,s`$ or $`b`$, respectively) or $`f=\tau `$ is the fermion in the loop and $`S=\stackrel{~}{\nu }_k`$ or $`\stackrel{~}{\nu }_k^{}`$ ($`k=1,2`$ or 3), is the scalar responsible for the flavor changing transition $`\mu e`$. Also

$`A_{RPV}^{\tau \stackrel{~}{\nu }_k}=\lambda _{k12}\lambda _{k33}^2,A_{RPV}^{\tau \stackrel{~}{\nu }_k^{}}=\lambda _{k21}\lambda _{k33}^2,A_{RPV}^{d_i\stackrel{~}{\nu }_k}=\lambda _{k12}\lambda _{kii}^{}^2/9,A_{RPV}^{d_i\stackrel{~}{\nu }_k^{}}=\lambda _{k21}\lambda _{kii}^{}^2/9.`$ (7)
The branching ratios scaled by the relevant RPV couplings for single particles running in the triangle are given in Table 1. In the spirit of this paper, only such down quarks and leptons running in the triangle are considered that would not give rise to the decay $`Ll\gamma `$, i.e., only those leptons that cannot generate the $`\lambda `$-reducible-type diagrams.

Comparing the current limit $`BR(\mu e\gamma \gamma )<7.2\times 10^{11}`$ PDG  with the estimate for the maximal possible branching ratio for this decay mode in Table 1, in conjunction with the existing bounds on the relevant RPV coupling products as given in Table 2, it can be argued that, at present, $`\mu e\gamma \gamma `$ does not impose any new constraints on the RPV SUSY parameter space. However, the situation will improve with the much more stringent experimental constraints that can be anticipated: on the third row of Table 2 the expected sensitivity of $`\mu e\gamma \gamma `$ to the relevant RPV coupling products is exhibited, assuming that future experiments will be sensitive to $`BR(\mu e\gamma \gamma )\stackrel{<}{}10^{14}`$, see e.g. experiment4 .

For the case of $`\tau `$-decays, taking the existing bounds on the relevant RPV coupling products into account, the largest branching ratio is obtained for $`\tau \mu \gamma \gamma `$ with an $`s`$-quark in the loop. In particular, for $`|\lambda _{233}\lambda _{322}^{}|=0.01`$ (its upper bound, see SK-02 ), $`m_s=450`$ MeV and $`M_{\stackrel{~}{\nu }}=100`$ GeV, $`BR(\tau \mu \gamma \gamma )10^{10}`$. This result is by far smaller than the typical present limits on $`\tau `$ decay branching ratios, which are of the order of $`10^7`$, see e.g. CLEO . Moreover, planned experiments Turks , which will produce about $`5\times 10^8`$ taus per year, will also be insensitive to the decays $`\tau l\gamma \gamma `$.

A few remarks are in order to conclude the discussion of the branching ratios:

* There is also a second possibility for $`Ll\gamma \gamma `$ of the $`\lambda `$-irreducible type to proceed: through the soft SUSY-breaking RPV terms ADL  $`C_{ijk}\stackrel{~}{L}_i\stackrel{~}{L}_j\stackrel{~}{E}_k^c`$ or $`C_{ijk}^{}\stackrel{~}{L}_i\stackrel{~}{Q}_j\stackrel{~}{D}_k^c`$ and slepton or squark in the loop. In this case the sneutrino couples to a pair of sleptons or down squarks of different helicity index, i.e. $`\stackrel{~}{\nu }\stackrel{~}{f}_L\stackrel{~}{f}_R`$, where $`\stackrel{~}{f}`$ is a slepton or down squark, while the photons couple to $`\stackrel{~}{f}_L\stackrel{~}{f}_L`$ or to $`\stackrel{~}{f}_R\stackrel{~}{f}_R`$. Assuming for simplicity that the $`\stackrel{~}{f}_L\stackrel{~}{f}_R`$ mixing is of $`𝒪(1)`$ and setting $`C_{ijk}M_{SUSY}\stackrel{~}{C}_{ijk}`$, $`C_{ijk}^{}M_{SUSY}\stackrel{~}{C^{}}_{ijk}`$, the width for $`Ll\gamma \gamma `$ is (for $`m_{\stackrel{~}{f}}^2/m_L^21`$ and $`m_l^2/m_L^21`$):

$`\mathrm{\Gamma }(Ll\gamma \gamma )\stackrel{~}{A}_{RPV}^{ij\stackrel{~}{f}S}\left({\displaystyle \frac{\alpha }{\pi }}\right)^2{\displaystyle \frac{1}{256\pi ^3}}{\displaystyle \frac{Q_{\stackrel{~}{f}}^4}{8640}}{\displaystyle \frac{m_L^7M_{SUSY}^2}{m_{\stackrel{~}{f}}^4m_S^4}},`$ (8)
where $`\stackrel{~}{f}=\stackrel{~}{\mathrm{}}`$ with $`Q_{\stackrel{~}{f}}=1`$ or $`\stackrel{~}{f}=\stackrel{~}{d}`$ with $`Q_{\stackrel{~}{f}}=1/3`$ if the sparticle in the loop is a slepton or down squark, respectively, and $`S=\stackrel{~}{\nu }`$ or $`\stackrel{~}{\nu }^{}`$. Also,

$`\stackrel{~}{A}_{RPV}^{ij\stackrel{~}{\mathrm{}}_m\stackrel{~}{\nu }_k}=\stackrel{~}{C}_{kmm}\lambda _{kij}^2,\stackrel{~}{A}_{RPV}^{ij\stackrel{~}{\mathrm{}}_m\stackrel{~}{\nu }_k^{}}=\stackrel{~}{C}_{kmm}\lambda _{kji}^2,`$
$`\stackrel{~}{A}_{RPV}^{ij\stackrel{~}{d}_m\stackrel{~}{\nu }_k}=9\stackrel{~}{C^{}}_{kmm}\lambda _{kij}^2,\stackrel{~}{A}_{RPV}^{ij\stackrel{~}{d}_m\stackrel{~}{\nu }_k^{}}=9\stackrel{~}{C^{}}_{kmm}\lambda _{kji}^2.`$ (9)
Here $`m,k=1,2,3`$ denote the family of the slepton ($`\stackrel{~}{\mathrm{}}`$) or down squark ($`\stackrel{~}{d}`$) and sneutrino or antisneutrino, respectively.

Plugging in numbers (e.g., setting $`M_{SUSY}=m_{\stackrel{~}{\mathrm{}}_m}=m_{\stackrel{~}{\nu }_k}=100`$ GeV), $`BR(\mu e\gamma \gamma )\stackrel{~}{C}_{kmm}\lambda _{k12}^23.7\times 10^{14}`$ and $`BR(\tau e\gamma \gamma )\stackrel{~}{C}_{kmm}\lambda _{k13}^21.9\times 10^{12}`$. Both values are much smaller than the ones obtained with fermions in the loop.
* If an RPV bilinear term exists, i.e. $`B\stackrel{~}{L}H_u`$ ADL , then $`Ll\gamma \gamma `$ can also proceed via the $`\lambda `$-irreducible topology with one trilinear ($`\lambda `$) and one bilinear ($`B`$) RPV insertion. Therefore, the decay would be $`\lambda \times B`$. In this case the off-shell sneutrino emitted from the decaying lepton mixes with the down-type Higgs boson, which then goes into the two photons with the well known (see e.g. HHG ) one-loop SUSY amplitude for $`H\gamma \gamma `$ (with the replacement $`m_H^22kk^{}`$). This possibility will be investigated in detail elsewhere AGt .
* Another set of related processes that may be generated through the same $`\lambda `$-irreducible diagrams are $`\tau (e,\mu )gg`$. Work on these processes is in progress AGt .
* The $`\tau `$-decays $`\tau (\mu ,e)\gamma \gamma `$ of the $`\lambda `$-irreducible type that are considered in this paper, with the lepton loop, have tree-level “analogs” that probe identical combinations of lambdas: $`\tau \mu ee`$ and $`\tau \mu \mu e`$ CR-96 . However, the contribution from lepton loops is not necessarily the largest one: it was found above that the dominant contribution stems from the $`s`$-quark loop, for which the branching ratio can amount to up to $`𝒪(10^{10})`$ (for $`\tau \mu \gamma \gamma `$), taking into account the current constraints SK-02 , which are obtained from processes of the type $`\tau l+\mathrm{hadrons}`$.
* As it can be seen from the Table, $`\mu e`$ conversion in nuclei is a strong rival to our processes. However, couplings that involve $`b`$ quark and $`\tau `$ lepton (and perhaps $`s`$ quark) are beyond the reach of $`\mu e`$ conversion experiments. Although our bounds on such processes are weak, they will strengthen if PRISM experiment will come into action.

In this paper a new class of lepton flavor-violating-processes, $`Ll\gamma \gamma `$ has been presented and discussed that might appear even if the corresponding decay into one photon only is suppressed or not existing at all. The underlying mechanism for such decays is the decay of the heavy lepton into the lighter one and a virtual scalar particle that goes into the two photons via a triangle diagram. Hence, such processes might probe the flavor structure of leptons coupling to these scalar particles that are a common feature of many models for new physics. Branching ratios for such decays have been calculated within the framework of R-parity-violating supersymmetry with sneutrinos as the scalar particles. Imposing RPV bounds emerging from other processes, the resulting branching ratios are well below current observation thresholds. However, for a planned new round of experiments, especially for $`\mu `$ decays, these processes - if they exist - might shed light on potential new physics or in turn might help to set more stringent exclusion bounds.

###### Acknowledgements.

A.G. acknowledges financial support by DAAD and thanks G. Soff and the members of his group at the Dresden Technical University for providing an atmosphere conducive to research during his visit in Dresden. S.B.S. would like to thank S. Roy for helpful discussions. This work is supported by the Israel Science Foundation and by the VPR Fund for the Promotion of Research at the Technion.


format: text

text:
# Light clusters production as a probe to nuclear symmetry energy

thanks: On leave from Department of Physics, Shanghai Jiao Tong University, Shanghai 200030, China

## Abstract

Using an isospin-dependent transport model for heavy-ion collisions induced by neutron-rich nuclei at intermediate energies, we study the production of light clusters such as deuteron, triton, and 3He via coalescence of nucleons. We find that both the yields and energy spectra of these light clusters are affected significantly by the density dependence of nuclear symmetry energy, with a stiffer symmetry energy giving a larger yield.

The energy per particle of an asymmetric nuclear matter with density $`\rho `$ and an isospin asymmetry $`\delta =(\rho _n\rho _p)/\rho `$, where $`\rho _n`$ and $`\rho _p`$ are, respectively, its neutron and proton densities, is usually approximated by a parabolic law li98 , i.e.,

$$E(\rho ,\delta )=E(\rho ,0)+E_{\mathrm{sym}}(\rho )\delta ^2,$$
(1)
where $`E(\rho ,0)`$ is the energy per particle of symmetric nuclear matter and $`E_{\mathrm{sym}}(\rho )`$ is the nuclear symmetry energy. Although the nuclear symmetry energy at normal nuclear matter density $`\rho _0=0.16`$ fm-3 has been determined to be around $`30`$ MeV from the empirical liquid-drop mass formula myers ; pomorski , its values at other densities are poorly known. Studies based on various theoretical models also give widely different predictions ibook . Lack of this knowledge has hampered our understanding of both the structure of radioactive nuclei oya ; brown ; hor01 ; furn02  and many important issues in nuclear astrophysics bethe ; bom ; lat01 , such as the nucleosynthesis during pre-supernova evolution of massive stars and the properties of neutron stars bethe ; lat01 . However, recent advance in radioactive nuclear beam facilities provides the opportunity to study the density dependence of the nuclear symmetry energy. Theoretical studies have already shown that in heavy-ion collisions induced by neutron-rich nuclei, the effect of nuclear symmetry energy can be studied via the pre-equilibrium neutron/proton ratio li97 , the isospin fractionation fra1 ; fra2 ; xu00 ; tan01 , the isoscaling in multifragmentation betty , the proton differential elliptic flow lis , the neutron-proton differential transverse flow li00 , the $`\pi ^{}`$ to $`\pi ^+`$ ratio li02 , and two-nucleon correlation functions chen02 .

In this work, we study the production of deuteron, triton, and 3He in heavy-ion collisions induced by neutron-rich nuclei by means of the coalescence model based on the nucleon phase space distribution functions from an isospin-dependent transport model. It is found that both the multiplicities and energy spectra of these light clusters are sensitive to the density dependence of nuclear symmetry energy but not to the isospin-independent part of nuclear equation of state and the in-medium nucleon-nucleon cross sections. Therefore, light clusters production in heavy-ion collisions induced by neutron-rich nuclei provides another possible method for extracting useful information about the nuclear symmetry energy.

The coalescence model has been used extensively in describing the production of light clusters in heavy-ion collisions at both intermediate Gyu83 ; Aich87 ; Koch90 ; Indra00  and high energies Mattie95 ; Nagle96 . In this model, the probability for producing a cluster is determined by its Wigner phase-space density and the nucleon phase-space distribution at freeze out. Explicitly, the multiplicity of a $`M`$-nucleon cluster in a heavy ion collision is given by Mattie95

$`N_M`$ $`=`$ $`G{\displaystyle 𝑑𝐫_{i_1}𝑑𝐪_{i_1}\mathrm{}𝑑𝐫_{i_{M1}}𝑑𝐪_{i_{M1}}}`$ (2)
$`\times `$ $`\underset{i_1>i_2>\mathrm{}>i_M}{{\displaystyle }}\rho _i^W(𝐫_{i_1},𝐪_{i_1}\mathrm{}𝐫_{i_{M1}},𝐪_{i_{M1}}).`$
In the above, $`𝐫_{i_1},\mathrm{},𝐫_{i_{M1}}`$ and $`𝐪_{i_1},\mathrm{},𝐪_{i_{M1}}`$ are, respectively, the $`M1`$ relative coordinates and momenta taken at equal time in the $`M`$-nucleon rest frame; $`\rho _i^W`$ is the Wigner phase-space density of the $`M`$-nucleon cluster; and $`\mathrm{}`$ denotes event averaging. The spin-isospin statistical factor for the cluster is given by $`G`$, and its value is $`3/8`$ for deuteron and $`1/3`$ for triton or 3He, with the latter including the possibility of coalescence of a deuteron with another nucleon to form a triton or 3He Polleri99 .

For the deuteron Wigner function, it is obtained from the Hulthén wave function, i.e.,

$$\varphi =\sqrt{\frac{\alpha \beta (\alpha +\beta )}{2\pi (\alpha \beta )^2}}\frac{e^{\alpha r}e^{\beta r}}{r}$$
(3)
with parameters $`\alpha =0.23`$ fm-1 and $`\beta =1.61`$ fm-1 to reproduce the measured deuteron root-mean-square radius of $`1.96`$ fm Ericsson88 . As in Ref. Mattie95 , we express the Hulthén wave function in terms of $`15`$ arbitrary Gaussian functions and determine their strengths and ranges by a least square fit. For triton and 3He, their Wigner functions are obtained from the product of three ground-state wave functions of a spherical harmonic oscillator with its parameters adjusted to reproduce the measured root-mean-square radii of triton and 3He, i.e., $`1.61`$ fm and $`1.74`$ fm, respectively chen86 . Normal Jacobian coordinates for a three-particle system are then introduced to derive the Wigner functions for triton and 3He as in Ref. Mattie95 .

The space-time distribution of nucleons at freeze out is obtained from an isospin-dependent Boltzmann-Uehling-Uhlenbeck (IBUU) transport model (e.g., li97 ; li00 ; li02 ; li96 ). For a review of the IBUU model, we refer the reader to Ref. li98 . The IBUU model includes explicitly the isospin degree of freedom through different proton and neutron initial distributions as well as their different mean-field potentials and two-body collisions in subsequent dynamic evolutions. For the nucleon-nucleon cross sections, we use as default the experimental values in free space. In order to study the effects due to the isospin-dependence of in-medium nucleon-nucleon cross sections $`\sigma _{\mathrm{medium}}`$, we also use a parameterization obtained from the Dirac-Brueckner approach based on the Bonn A potential lgq9394 . For experimental free-space cross sections, the neutron-proton cross section is about a factor of $`3`$ larger than the neutron-neutron or proton-proton cross sections. On the other hand, the in-medium nucleon-nucleon cross sections used here have smaller magnitudes and weaker isospin dependence than $`\sigma _{\mathrm{exp}}`$ but strong density dependence. For the isoscalar potential, we use as default the Skyrme potential with an incompressibility $`K_0=380`$ MeV. This potential has been shown to reproduce the transverse flow data from heavy ion collisions equally well as a momentum-dependent soft potential with $`K_0=210`$ MeV pan93 ; zhang93 .

The IBUU model is solved with the test particle method buu . Although the mean-field potential is evaluated with test particles, only collisions among nucleons in each event are allowed. Light cluster production from coalescence of nucleons is treated similar as nucleon-nucleon collisions, i.e., only nucleons in the same event are allowed to coalesce to light clusters. Results presented in the following are obtained with $`20,000`$ events using $`100`$ test particles for a physical nucleon.

For the density dependence of symmetry energy, we adopt the parameterization used in Ref. hei00  for studying the properties of neutron stars, i.e.,

$$E_{\mathrm{sym}}(\rho )=E_{\mathrm{sym}}(\rho _0)u^\gamma ,$$
(4)
where $`u\rho /\rho _0`$ is the reduced density and $`E_{\mathrm{sym}}(\rho _0)=30`$ MeV is the symmetry energy at normal nuclear matter density. In the following, we consider the two cases of $`\gamma =0.5`$ (soft) and $`2`$ (stiff) to explore the large range of $`E_{\mathrm{sym}}(\rho )`$ predicted by many-body theories bom .

We consider the reaction of 52Ca + 48Ca, which has an isospin asymmetry $`\delta =0.2`$ and can be studied at future Rare Isotope Accelerator Facility. In the present study, nucleons are considered as being frozen out when their local densities are less than $`\rho _0/8`$ and subsequent interactions do not cause their recapture into regions of higher density. Shown in Figs. 1 (a), (b) and (c) are time evolutions of the average multiplicity of deuteron, triton, and 3He from central collisions of 52Ca + 48Ca at $`E=80`$ MeV/nucleon by using the soft (solid curve) and stiff (dashed curve) symmetry energies. It is seen that production of these light clusters from a neutron-rich reaction system is very sensitive to the density dependence of nuclear symmetry energy. Final multiplicities of deuteron, triton, and 3He for the stiff symmetry energy is larger than those for the soft symmetry energy by $`51\%`$, $`73\%`$, and $`100\%`$, respectively. This is due to the fact that the stiff symmetry energy induces a stronger pressure in the reaction system and thus causes an earlier emission of neutrons and protons than in the case of the soft symmetry energy chen03 , leading to a stronger correlations among nucleons. Furthermore, the soft symmetry energy, which gives a more repulsive symmetry potential for neutrons and more attractive one for protons in low density region ($`\rho _0`$) than those from the stiff symmetry energy, generates a larger phase-space separation between neutrons and protons at freeze out, and thus a weaker correlations among nucleons. The larger sensitivity of the multiplicity of 3He to the nuclear symmetry energy than those of triton and deuteron as seen in Fig. 1 reflects the fact that the symmetry energy effect is stronger on lower momentum protons than neutrons chen03 .

Although the density dependence of nuclear symmetry energy affects appreciably the yield of light clusters, changing the incompressibility from $`K_0=380`$ to $`201`$ MeV (dotted curves) or using $`\sigma _{\mathrm{medium}}`$ instead of $`\sigma _{\mathrm{exp}}`$ (dash-dotted curves) only leads to a small change in the yield of these clusters as shown in Fig. 1. This implies that the relative space-time structure of neutrons and protons at freeze out is not sensitive to the equation of state (EOS) of symmetric nuclear matter and in-medium nucleon-nucleon cross sections.

The kinetic energy spectra in the center-of-mass system for deuteron, triton, and 3He are shown in Fig. 2 for both the soft (solid squares) and stiff (open squares) symmetry energies. It is seen that the symmetry energy has a stronger effect on the yield of low energy clusters than that of high energy ones. For example, the symmetry energy effect on the yield of deuteron, triton, and 3He is about $`60\%`$, $`100\%`$, and $`120\%`$, respectively, if their kinetic energies are around $`10`$ MeV, but is about $`30\%`$, $`40\%`$, and $`85\%`$, respectively, if their kinetic energies are around $`100`$ MeV. This follows from the fact that lower energy clusters are emitted later in time when the size of nucleon emission source is relatively independent of nuclear symmetry energy, leading thus to a similar probability for nucleons to form light clusters. Since there are more low energy nucleons for a stiffer symmetry energy, more light clusters are thus produced. On the other hand, higher energy nucleons are emitted earlier when the size of emission source is more sensitive to the symmetry energy, with a smaller size for a stiffer symmetry energy. The probability for light cluster formation is thus larger for a stiffer symmetry energy. This effect is, however, reduced by the smaller number of high energy nucleons if the symmetry energy is stiffer. As a result, production of high energy light clusters is less sensitive to the stiffness of symmetry energy. This is different from that seen in the correlation functions between two nucleons with low relative momentum, where the symmetry energy effect is larger for nucleon pairs with higher kinetic energies chen02 , as they are only affected by the size of the emission source, not the number of emitted nucleon pairs.

The isobaric yield ratio t/3He is less model-dependent and also less affected by other effects, such as the feedback from heavy fragment evaporation and the feed-down from produced excited triton and 3He states. In Fig. 3, we show the t/3He ratio with statistical errors as a function of cluster kinetic energy in the center-of-mass system for the soft (solid squares) and stiff (open squares) symmetry energies. It is seen that the ratio t/3He obtained with different symmetry energies exhibits very different energy dependence. While the t/3He ratio increases with kinetic energy for the soft symmetry energy, it decreases with kinetic energy for the stiff symmetry energy. The symmetry energy thus affects more strongly the ratio of high kinetic energy triton and 3He. For both soft and stiff symmetry energies, the ratio t/3He is larger than the neutron to proton ratio of the whole reaction system, i.e., N/Z$`=1.5`$. This is in agreement with results from both experiments and statistical model simulations for other reaction systems and incident energies Cibor00 ; Hagel00 ; Vesel01 ; Chomaz99 . It is interesting to note that although the yield of lower energy triton and 3He is more sensitive to symmetry energy than higher energy ones, as shown in Fig. 2, their ratio at higher energy is affected more by the symmetry energy. Moreover, the energy dependence of the ratio t/3He is insensitive to the EOS of symmetric nuclear matter and in-medium nucleon-nucleon cross sections. These features thus imply that the pre-equilibrium triton to 3He ratio is also a sensitive probe to the density dependence of nuclear symmetry energy.

Isospin effects on cluster production and isotopic ratios in heavy ion collisions have been previously studied using either the lattice gas model Chomaz99  or a hybrid of IBUU and statistical fragmentation model tan01 . These studies are, however, at lower energies than considered here, where effects due to multifragmentation as a result of possible gas-liquid phase transition may play an important role. Except deuterons, both tritons and 3He are only about one per event in heavy ion collisions at $`80`$ MeV/nucleon. The number of other clusters such as the alpha particle is not large either jacak ; borderie . In this case, the coalescence model is expected to be a reasonable model for determining the production of light clusters from heavy ion collisions. Furthermore, the effect obtained in present study will be enhanced if other clusters are emitted earlier as the isospin asymmetry of the residue is increased.

In conclusion, using an isospin-dependent transport model together with a coalescence model for light cluster production, we have found that the nuclear symmetry energy affects significantly the production of light clusters in heavy-ion collisions induced by neutron-rich nuclei. More deuterons, tritons, and 3He are produced with the stiff nuclear symmetry energy than the soft nuclear symmetry energy. This effect is particularly large when these clusters have lower kinetic energies. Also, the isobaric ratio t/3He, especially for higher energy tritons and 3He, shows a strong sensitivity to the density dependence of nuclear symmetry energy. It is further found that light clusters production is not sensitive to the isospin-independent part of nuclear equation of state and the in-medium nucleon-nucleon cross sections. The study of light clusters production in heavy ion collisions induced by neutron-rich nuclei thus allows us to extract useful information about the density-dependence of nuclear symmetry energy.

We thank Joe Natowitz for discussions on light clusters production in heavy ion collisions. This paper is based on the work supported by the U.S. National Science Foundation under Grant Nos. PHY-0098805 and PHY-0088934 as well as the Welch Foundation under Grant No. A-1358. LWC is also supported by the National Natural Science Foundation of China under Grant No. 10105008.


text:
# Evolutionary models for cool brown dwarfs and extrasolar giant planets. The case of HD 209458

## 1 Introduction

The past decade was marked by two major discoveries in the field of stellar and planetary physics: the detections of the first unambiguous brown dwarf (BD) GL 229B (Oppenheimer et al . 1995) and the first extrasolar giant planet (EGP) 51 Peg b (Mayor & Queloz 1995). The near-IR spectrum of GL 229B was found to be dominated by strong methane absorption bands, looking more similar to Jupiter than to late type-stars. On the other hand, the surprisingly small orbital separation between 51 Peg b and its parent star suggests that the planet should be affected by irradiation and that, given the expected large surface temperature, its atmospheric properties should resemble more the ones of relatively hot brown dwarfs than the ones of jovian planets.

Since then, about thirty methane dwarfs (or the so-called T-dwarfs) have been identified, due mainly to the near-IR surveys 2MASS (Burgasser et al. 1999), SDSS (Strauss et al. 1999) and the VLT (Cuby et al. 1999). The radial velocity technique has now revealed more than 100 EGPs in orbit around nearby stars (see Hubbard, Burrows, Lunine, 2002 for a review and references therein), with a large fraction ($``$ 10%-20%) being extremely close (less than 0.06 AU) to their parent star. The mass of substellar companions detected by radial velocimetry extends well above the deuterium burning minimum mass $`0.012M_{}=12M_\mathrm{J}`$ (Saumon et al., 1996; Chabrier et al., 2000a). This mass is often used as the boundary between planets and brown dwarfs, more for semantical than physical reasons. On the other hand isolated objects with planetary masses are now discovered in young stellar clusters, down to a few $`M_\mathrm{J}`$, as recently reported by Zapatero et al. (2002) in $`\sigma `$ Orionis. These observations suggest that there is an overlap between the mass range of the least massive brown dwarfs and of the most massive giant planets. In principle, different formation processes should distinguish planets from brown dwarfs. However, such a distinction is difficult to characterize in terms of atmospheric, structural and cooling properties since both types of objects have convective interiors with essentially a metallic H/He mixture. The signature of a central rock+ice core, like in solar giant planets, would be the clear identification of a planet. The presence of a core can affect the radius of a planet, yielding a smaller planetary radius than in the absence of a core. For 1 $`M_\mathrm{J}`$, the effect is about 5% on the radius for a core mass $`<`$ 0.06 $`M_\mathrm{J}`$ (see Saumon et al. 1996). The presence of this core can be inferred from the accurate characterization of the gravitational moments of the object, and such an observation is currently not feasible for EGPs. In addition, both giant planets and brown dwarfs have atmospheres dominated by molecular absorption and the effects of cloud formation. Although frustrating from an observational point of view, these similarities imply that the general cooling theory for BDs, involving detailed models of the atmosphere and inner structures, can be applied to EGPs. In terms of cooling properties, this general theory can even be applied to Jupiter, as emphasized by Hubbard et al. (2002). Additional observational constraints, as provided by spacecraft encounters or by direct probes (e.g Galileo), have lead to refinements of the models (heavy element core, non-standard chemical composition). As mentioned above, such constraints are, unfortunately, far from being accessible for EGPs.

Much effort has been devoted to the modeling of substellar objects during the past decade, improving our understanding of cool atmospheres (see Allard et al. 1997 for a review), of the role of dust (Tsuji et al. 1999; Burrows et al. 2000; Ackerman & Marley 2001, Allard et al. 2001; Marley et al. 2002), of irradiation (Saumon et al. 1996; Seager & Sasselov 1998; Sudarsky et al. 2000; Barman et al. 2001), and of their inner structure and evolutionary properties (Burrows et al. 1997 ; Chabrier et al. 2000b; see Chabrier & Baraffe 2000 for a review). One remaining major challenge in the theory is the description of dynamical processes of grain formation and diffusion necessary to understand the transition between L-dwarfs and T-dwarfs, which is expected to take place at $`T_{\mathrm{eff}}`$ 1300K - 1700K. The former objects are better reproduced by dusty atmosphere models, whereas the later are better reproduced by dust-free (or partly dusty) models. The recent observations of L/T dwarfs at the transition clearly indicate that complex processes take place in the atmosphere of these objects (see e.g Burgasser et al. 2002). Another important challenge is the modeling of irradiation effects, which are expected to affect the spectra of close-in EGPs, and may also affect their inner structure and cooling properties. The recent discovery of the planet HD209458b transiting its parent star (Charbonneau et al. 2000) provides a unique test to explore such effects, since its mass and radius can be determined with high accuracy from the modeling of the transit lightcurve. According to the most recent determination (Cody & Sasselov 2002), the mass and radius of the planet are estimated to be $`m=0.69\pm 0.02M_\mathrm{J}`$ and $`R=1.42_{0.13}^{+0.10}R_\mathrm{J}`$.

Evolutionary models including crude estimates of the effects of irradiation on planet atmospheres suggest that extrinsic heating is sufficient to maintain a larger planetary radius compared to an isolated planet. It has thus been argued that irradiation could explain the large radius of HD209458b (Guillot et al. 1996; Burrows et al. 2000). More recently, Guillot & Showman (2002) questioned such results and argue that the radius of HD209458b can only be reproduced if the deep atmosphere is much hotter than what can be expected from irradiation effects. However, none of these calculations includes a consistent treatment between the irradiated atmospheric structure and the interior structure of the planet. Such a consistent treatment is mandatory to get reliable results since the deep interior entropy profile, which determines the heat content of the planet to be radiated away while it cools, is affected by the modification of the atmospheric temperature profile due to the incoming external heat flux. The main goal of the present paper is to present the first such consistent calculations. As mentioned above, in the case of non-irradiation, these calculations apply to the evolution of cool (dust-free like) brown dwarfs, i.e. T-dwarfs, and extrasolar giant planets far enough from the parent star for the irradiation effects on the thermal structure of the planet to be negligible. This is the case of the solar giant planets, the cooling of which is simply characterized by the cooling properties of the ”isolated” planet plus the heating contribution from the Sun $`4\pi \sigma R_p^2T_{}^4`$, where $`R_p`$ is the radius of the planet and $`T_{}`$ represents the equivalent black body temperature of the converted solar radiation (Hubbard, 1977; Guillot et al., 1995). In section 2 we briefly present the input physics of non-irradiated models, describing methane dwarfs and isolated EGPs. Apart from the impinging stellar flux, the same input physics are used to analyse the effects of irradiation. The effects on the radius and cooling properties of giant planets are described in section 3 and results are compared to the observed properties of HD209458b. Discussion follows in §4.

## 2 Non-irradiated models

### 2.1 Model description

The main input physics involved in the present models are the same as described in our previous works (Chabrier & Baraffe 1997; Baraffe et al. 1998; Chabrier et al. 2000b). The models are based on the coupling between interior and non-grey atmosphere structures. The treatment of dust in the atmosphere is described in detail in Allard et al. (2001), with two limiting cases of dust treatment. The first case, referred to as “DUSTY”, takes into account the formation of dust in the equation of state, and its scattering and absorption in the radiative transfer equation. Such models assume that dust species remain where they form, according to the chemical equilibrium conditions. The second case, referred to as “COND”, neglects dust opacity in the radiative transfer equation. In a previous paper (Chabrier et al. 2000b), we presented the evolutionary models based on DUSTY atmosphere models, aimed at describing the evolution and the photometric and spectroscopic properties of early L-dwarfs. The present paper is devoted to evolutionary models based on the COND approach, which are more appropriate to objects with effective temperatures $`T_{\mathrm{eff}}<\mathrm{\hspace{0.17em}1300}`$K, such as methane dwarfs or EGPs at large orbital separation. These models apply when all grains have gravitationally settled below the photosphere. A forthcoming paper will be devoted to models taking into account characteristic diffusion timescales of different processes affecting the dust stratification (e.g coagulation, gravitational settling, convection). These models aim at describing in particular the transition objects between late L-dwarfs and early T-dwarfs (see Allard 2002).

### 2.2 Evolution of cool brown dwarfs

A preliminary version of the COND models was presented in Chabrier et al. (2000b) down to 0.01 $`M_{}`$. In the present paper, we extend our calculations to $`T_{\mathrm{eff}}`$ = 100 K and $`m=`$ 0.5 $`M_\mathrm{J}`$. The evolution of $`L`$ and $`T_{\mathrm{eff}}`$ as a function of time for different masses is displayed in Fig. 1. The properties of the COND models for different ages are given in Tables 1-5. As already stressed in Chabrier et al. (2000b), the treatment of dust in the atmospheric models barely affects the evolutionary tracks in terms of $`L`$ and $`T_{\mathrm{eff}}`$ as a function of time for a given mass (see Fig. 2 of Chabrier et al. 2000b). Consequently, although the COND models are more appropriate to describe the spectral and photometrical properties of substellar objects with $`T_{\mathrm{eff}}<`$ 1300K, they provide a good description of their cooling properties even at higher $`T_{\mathrm{eff}}`$. In other words, it is not necessary to perform evolutionary calculations with the DUSTY models above $`T_{\mathrm{eff}}`$ 1300K and switch to the COND models below. An object characteristic of the present models was recently discovered by Zapatero et al. (2002) in $`\sigma `$ Orionis: S Ori 70. From a comparison of its observed spectrum with COND synthetic spectra (Allard et al. 2001), Zapatero et al. (2002) estimate an effective temperature $`T_{\mathrm{eff}}`$ 700-1000K. If the membership of S Ori 70 to $`\sigma `$ Orionis is confirmed, implying an age $`<`$ 10 Myr, its mass mass should be $`<\mathrm{\hspace{0.17em}5}M_\mathrm{J}`$ (see Fig. 1).

Several methane dwarfs have been discovered in the solar neighbourhood, implying older ages and thus larger masses than the extreme case of S Ori 70. At an age of $`10^8`$ yr, only objects with masses below the deuterium burning minimum mass ($`m0.012M_{}`$) have $`T_{\mathrm{eff}}<\mathrm{\hspace{0.17em}1300}K`$, whereas at 5 Gyr, it is the case for all substellar objects with $`m<\mathrm{\hspace{0.17em}0.06}M_{}`$ (see Fig. 1 and Tables 1 and 4). Photometric observations and parallax determinations of several L- and T-dwarfs (Els et al. 2001; Leggett et al. 2002a; Dahn et al. 2002) now allow a comparison with models in observational color-magnitude diagrams (CMD), providing stringent constraints on theoretical models (see Figs. 2-4).

As already noticed in Allard et al. (1996) for GL 229B, models free of atmospheric dust clouds better reproduce the near-IR photometric and spectral properties of methane dwarfs. This is illustrated in $`(JK)`$ and $`(KL^{})`$ colors in Fig. 2 and Fig. 3 respectively, where the COND models reproduce the main observed trends. In Fig. 2, we note the two transition objects, intermediate between L- and T- dwarfs, with $`(JK)1`$ (GL 86B: Els et al. 1999; SDSS 1254-01: Leggett et al. 2002a, Dahn et al. 2002) and the faintest L-dwarfs (Dahn et al. 2002), which are not described by either the DUSTY or COND limiting cases, and require a detailed treatment of dust diffusion in the atmosphere, as mentioned in §2.1.

The predictions of the COND models provide a general good agreement with observed near-IR photometry and spectra at wavelength $`>1\mu `$m (Leggett et al. 2002b; Zapatero et al. 2002). The models show however shortcomings at shorter wavelength, with a flux excess around 0.8-0.9 $`\mu `$m, characteristic of the I-bandpass. This problem is illustrated in Fig. 4 in a $`(IJ)`$ \- $`M_\mathrm{J}`$ CMD, where the COND models predict significantly bluer $`(IJ)`$ colors than observations. As mentioned in Allard et al. (2001), uncertainties in the current treatment of the far wings of the absorption lines of alkali elements (Na, K) at such pressures may be responsible for this discrepancy. No theory, however, exists to date for an accurate description of broadening of atomic lines by collisions with H2 and He. Attempts to improve current treatments are under progress (Burrows & Volobuyev 2002).

The correct trend of colors and spectral properties predicted by the present models at wavelength $`>1\mu `$m, where most of the flux is emitted for the concerned range of $`T_{\mathrm{eff}}`$, comfort us however with their reliability to describe extremely cool objects.

## 3 Irradiated models

### 3.1 Effect on atmosphere structure

As mentioned in the introduction, a non-negligible fraction of EGPs orbit close to their parent star and their thermal and mechanical structure is affected by irradiation effects. Therefore, a general theory of cool substellar objects must take these effects into account. Recently Barman et al. (2001) have modeled irradiated atmospheres by including the impinging radiation field in the solution of the radiative transfer equation. As shown by these authors, for a given intrinsic luminosity, non-irradiated planets have very different temperature structures than irradiated planets. Thus, substituting non-irradiated atmospheric structures with $`T_{\mathrm{eff}}`$ = $`T_{eq}`$ (see definition below, Eq.(8)) for irradiated structures, as done up to now in the literature, yields incorrect inner boundary conditions for evolutionary calculations (see e.g. figure 13 of Barman et al. (2001)). Given the present lack of an accurate treatment of atmospheric dust diffusion, the calculations were performed only for the DUSTY and COND limit cases, respectively. The results emphasize the strong dependence of the emergent spectrum and atmospheric structure on the presence or absence of dust. In the absence of dust, the impinging flux can penetrate in deeper layers of the planet atmosphere, affecting more drastically the inner structure of the planet than in the dusty case.

Except for a possible detection of sodium absorption in the atmosphere of HD 209458b (Charbonneau et al. 2002), no constraints on the atmospheric composition of EGPs are available at the present time. The only strong observational constraint available for irradiated models is the transit planet HD 209458b. The determination of its mass and radius provides a stringent test to irradiated atmosphere calculations and to the resulting structure and evolution. We thus apply our calculations of irradiated EGPs to HD 209458-like systems.

We have computed a grid of irradiated atmosphere models based on the COND input physics described in §2, as in Barman et al. (2001). Although more appropriate for EGPs with $`T_{\mathrm{eff}}<\mathrm{\hspace{0.17em}1300}`$ K, the COND models maximise the effect of irradiation on the inner atmosphere structure and thus on the evolution of EPGs (Barman et al. 2001). The grid covers a wide range of $`T_{\mathrm{eff}}`$ from 40K to 100K, in steps of 20K, and from 100K to 2800K, in steps of 100K. It covers a range of surface gravities from $`\mathrm{log}g=2.5`$ to $`\mathrm{log}g=4.5`$, in steps of 0.5 dex. We adopt the characteristics of HD 209458, assuming for the primary an effective temperature $`T_{\mathrm{eff}}`$ = 6000 K, a radius $`R_{}`$ = 1.18 $`R_{}`$ (Mazeh et al. 2000; Cody & Sasselov 2002) and an orbital separation $`a`$ = 0.046 AU (Charbonneau et al. 2000). As in Barman et al. (2001), we make the simplifying assumptions that the impinging radiation field is isotropic and the incident flux $`F_{\mathrm{inc}}`$ is redistributed only over the dayside, i.e.

$$F_{\mathrm{inc}}=\frac{1}{2}(\frac{R_{}}{a})^2F_{},$$
(1)
where $`F_{}`$ is the total flux from the primary (see discussion in §4).

Before proceeding any further, we briefly re-specify definitions of fluxes (see e.g. Brett & Smith 1993), since use of various terminologies leads to confusion. In all cases, the integrated net flux $`F_{\mathrm{net}}`$, obtained from the solution of the transfer equation, is the intensity integrated over both in-coming and out-going directions $`\mu `$ ($`\mu =\mathrm{cos}\theta `$, where $`\theta `$ is the angle of incidence). Assuming there is no extra source or sink of energy (e.g no horizontal energy transfer), energy conservation implies that all the incident energy coming in must go out. Therefore, in the case of irradiation, the in-coming flux from the parent star cancels out the extra out-going, absorbed and reradiated flux due to the heating of the upper layers of the planet atmosphere (see Fig. 5). The in-coming flux at the surface is $`F_{\mathrm{in}}`$ = - $`F_{\mathrm{inc}}`$ and the out-going flux at the surface is $`F_{\mathrm{out}}`$ = $`F_{\mathrm{inc}}`$ \+ $`\sigma T_{\mathrm{eff}}^4`$, where $`\sigma T_{\mathrm{eff}}^4`$ defines the intrinsic, unperturbed flux $`\sigma T_{\mathrm{eff}}^4`$ of the initial, non-irradiated atmosphere structure. Energy conservation thus implies:

$$F_{\mathrm{net}}=F_{\mathrm{out}}+F_{\mathrm{in}}=\sigma T_{\mathrm{eff}}^4$$
(2)
The non-irradiation case ($`F_{\mathrm{inc}}=0`$) corresponds to the usual condition $`F_{\mathrm{net}}=F_{\mathrm{out}}=\sigma T_{\mathrm{eff}}^4`$.

Our atmosphere models, irradiated or not, are thus characterized by the parameters $`T_{\mathrm{eff}}`$ and $`g`$. Of course, the same net flux $`F_{\mathrm{net}}`$ corresponds to two different atmospheric structures, in the non-irradiated and irradiated case, because of the extra energy source $`F_{\mathrm{in}}0`$ in the latter case (see Fig. 5). Given the above definitions, the net flux characterizes the intrinsic luminosity, i.e. the rate of energy released by the planet as it contracts and cools down:

$$L_{int}=4\pi R_\mathrm{p}^2\sigma T_{\mathrm{eff}}^4=T\frac{dS}{dt}dm.$$
(3)
This quantity determines the cooling properties of the planet for a given set of outer boundary conditions provided by the atmospheric profile (see §3.2 below). We stress that, in the case of irradiation, $`T_{\mathrm{eff}}`$ does not characterize the total flux emitted by the planet, which is given by:

$$F_{\mathrm{out}}=\sigma T_{\mathrm{eff}}^4+F_{\mathrm{inc}}=\sigma T_{\mathrm{eff}}^4+\frac{1}{2}(\frac{R_{}}{a})^2F_{}$$
(4)
Note that $`F_{\mathrm{out}}`$ is the important quantity for observers, since it characterises the total radiation of the planet, including both thermal and reflected parts of the flux. However, we do not focus on this quantity, since a forthcoming paper will be devoted to spectral properties of irradiated planets (Barman et al. 2003, in preparation).

For the sake of comparison with non-irradiated atmosphere profiles (see Barman et al. 2001), we also define the quantity $`T_{\mathrm{therm}}`$ which characterises the thermal flux $`\sigma T_{\mathrm{therm}}^4`$ emitted by the irradiated fraction of the planet (in the present case, the day side only). This quantity reads:

$$\sigma T_{\mathrm{therm}}^4=\sigma T_{\mathrm{eff}}^4+(1A)F_{\mathrm{inc}},$$
(5)
where $`A`$ is the Bond albedo. According to the definitions above:

$$F_{\mathrm{out}}=\sigma T_{\mathrm{therm}}^4+F_{\mathrm{refl}},$$
(6)
where

$$F_{\mathrm{refl}}=AF_{\mathrm{inc}}$$
(7)
is the reflected part of the incident flux.

Within the conditions of the present calculations (Eq(1), $`a=0.046`$, $`T_{\mathrm{eff}}`$ = 6000 K, $`R_{}`$ = 1.18 $`R_{}`$), our Bond albedo is close to 0.1 for the coolest models ($`T_{\mathrm{eff}}100`$K). A final quantity, often used in the literature, is the equilibrium temperature, $`T_{\mathrm{eq}}`$, which characterizes the planet’s luminosity after having exhausted all its internal heat content (see e.g Guillot et al. 1996; Saumon et al. 1996):

$$\begin{array}{c}T_{\mathrm{eq}}^4=\frac{(1A)}{\sigma }F_{\mathrm{inc}}=\frac{(1A)}{2\sigma }(\frac{R_{}}{a})^2F_{}\hfill \\ \\ T_{\mathrm{therm}}^4\mathrm{when}T_{\mathrm{eff}}0.\hfill \end{array}$$
(8)
Note that given our definition of $`F_{\mathrm{inc}}`$, $`T_{\mathrm{eq}}^4`$ defined by Eq. (8) differs by a factor 2 from the definition usually used in the literature, because of the redistribution only over the day side. Note also that $`T_{\mathrm{eq}}`$ and $`T_{\mathrm{therm}}`$ differ significantly at young ages, when the intrinsic flux of the planet is not negligible.

The effect of irradiation on atmosphere structures is illustrated in Fig. 5 for different values of the effective temperature $`T_{\mathrm{eff}}`$. As already stressed in Barman et al. (2001), an irradiated structure characterised by $`T_{\mathrm{therm}}`$ can differ significantly from a non-irradiated structure at the same effective temperature $`T_{\mathrm{eff}}`$ = $`T_{\mathrm{therm}}`$. This point (see also Seager & Sasselov 1998; Guillot & Showman 2002) emphasizes the fact that adopting outer boundary conditions, for evolutionary calculations, from atmospheric profiles of nonirradiated models with $`T_{\mathrm{eff}}`$ = $`T_{\mathrm{therm}}`$, or $`T_{\mathrm{eff}}`$ = $`T_{\mathrm{eq}}`$ (as e.g. Burrows et al. (2000)), is incorrect and yields erroneous evolutionary properties for irradiated objects.

### 3.2 Effect on evolution

The main effect of irradiation on convective atmospheres and its consequences on evolution is well known (see Hubbard 1977; Brett & Smith 1993; Guillot et al 1996; Hubbard et al. 2002). The heating of the outer layers by the incident flux reduces the temperature gradient between these layers and the interior. They become radiative and the top of the convective zone is displaced to larger depths compared to the non-irradiated case, as clearly illustrated in Fig. 5. The inner atmosphere structure is hotter at a given pressure than the nonirradiated atmosphere model of same $`T_{\mathrm{eff}}`$ (see Fig. 5). In order to match the same inner entropy, or the same values of $`P`$ and $`T`$, characteristic of the boundary layer between the interior structure and the irradiated atmosphere structure, characterized by a given $`T_{\mathrm{eff}}`$ and $`\mathrm{log}g`$, one would need a nonirradiated atmosphere model with higher $`T_{\mathrm{eff}}`$, i.e a larger heat loss. Therefore, for a given planet heat content, i.e. internal entropy, the heat loss is reduced in the case of irradiation and the planet maintains a higher entropy for a longer time. Since for a given mass, the interior (P,T) profile and thus the entropy fix the radius, the irradiated planet has a larger radius than the nonirradiated counterpart at a given time, starting from the same initial configuration. In other terms, gravitational contraction, which is the dominant source of energy of the planet, proceeds more slowly with irradiation than without it.

Our calculations proceed as for our low-mass star or brown dwarf calculations, by coupling the interior and atmosphere profile at a deep enough optical depth, which defines unequivocally the fundamental properties of the object, $`m,R,T_{\mathrm{eff}},L`$ along its evolution $`t`$ (Chabrier & Baraffe, 1997). The boundary condition between inner and atmosphere structure is fixed at $`\tau _{\mathrm{std}}=100`$, which corresponds to a range of pressure $`P`$ = 0.1 - 200 bar for the whole range of atmosphere models 111$`\tau _{\mathrm{std}}`$ is defined at $`\lambda =1.2\mu m`$. The irradiated atmosphere models are integrated down to an optical depth $`\tau _{\mathrm{std}}=100`$ for $`T_{\mathrm{eff}}1000K`$ and $`\tau _{\mathrm{std}}=10^5`$ for $`T_{\mathrm{eff}}<`$ 1000K. In both cases, this is deep enough to reach the top of the convective zone and to provide a good spatial resolution of these layers, even for the coolest models (see Fig. 5). In any case, the incident flux $`F_{\mathrm{inc}}`$ drops to zero at $`\tau _{\mathrm{std}}<<50`$, well above the deepest layers of the atmosphere models. Note that for the coolest atmospheric structures, convection does not reach the layers corresponding to $`\tau _{\mathrm{std}}=100`$ (see Fig. 5). In that case, the radiative gradient in the interior is calculated with the Rosseland means of the same atmospheric opacities, for a consistent treatment between the interior and atmosphere thermal structures.

We have calculated the evolution of planets covering a range of mass from 0.5 $`M_\mathrm{J}`$ to 10 $`M_\mathrm{J}`$ with and without irradiation. The evolution of the radius as a function of time is shown in Fig. 6 for irradiated and nonirradiated EGPs of 1 $`M_\mathrm{J}`$ and 10 $`M_\mathrm{J}`$. As expected, the less massive the planet, the larger the effect of irradiation, for a given incident flux. At 1 Gyr, the 0.5 $`M_\mathrm{J}`$ EGP has a 14% larger radius than its nonirradiated counterpart, whereas for the 1 $`M_\mathrm{J}`$ (resp. 10 $`M_\mathrm{J}`$), $`R`$ is only 10% larger (resp. 7%). We also compare our COND models (the nonirradiated models) to the Burrows et al. (1997, hereafter B97) nonirradiated models. Significant differences appear at young ages ($`<1`$ Gyr), due certainly to different initial conditions (see Baraffe et al. 2002). For ages $`>1`$ Gyr and $`m>\mathrm{\hspace{0.17em}5}M_\mathrm{J}`$, the differences between the B97 models and ours are of the order of the irradiation effects. This reflects the different input physics, mainly in the dust treatment and molecular opacities and illustrates the present uncertainties in the models. For $`m<\mathrm{\hspace{0.17em}5}M_\mathrm{J}`$, however, irradiation effects become larger than the differences between the B97 and our models.

The specific case of HD209458b, with a mass $`m_\mathrm{p}=0.69M_\mathrm{J}`$, is illustrated in Fig. 7. The intrinsic luminosity and corresponding effective temperature in the irradiated case (solid lines) are compared to the non-irradiated case (dashed lines). Starting from the same initial configuration in both cases, the heat loss is reduced at early ages in the case of irradiation, as expected. Consequently, the irradiated model evolves at larger entropy and radius than its non-irradiated counterpart. During the first Myr of evolution, both evolutionary sequences contract with increasing central density and temperature, the non-irradiated model being denser. The latter becomes partially degenerate earlier, its contraction slows down and its heat loss becomes smaller than in the irradiated case (at $`\mathrm{log}t`$ 6.2 yr). The situation reverses at $`\mathrm{log}t`$ 7.4 yr when the effect of partial degeneracy becomes important in the irradiated sequence. The age of HD209458 is about 4-7 Gyr, according to Cody & Sasselov (2002). At 5 Gyr, the irradiated sequence displayed in Fig. 7 predicts a radius $`R=1.09R_\mathrm{J}`$, 26% smaller than the observed value. Without including irradiation effects, the radius is $`>`$ 30% than the observationally determined one. Note that the nonirradiated sequence stops at $`T_{\mathrm{eff}}`$ = 100K, corresponding to an age of $``$ 2 Gyr and a radius $`R=1R_\mathrm{J}`$. In the following section, we analyse the possible reasons for such a discrepancy.

## 4 Discussion

### 4.1 Uncertainties of irradiated atmosphere/evolutionary models

The question rises whether uncertainties of current models can explain the mismatch of HD209458b predicted versus observed radius, and whether irradiation effects can still provide the solution to the problem. We first note that our choice of parameters for the irradiated atmosphere calculations certainly overestimates the effects of irradiation (see §3). We assume redistribution of the incident flux over the day-side of the planet only, without taking into account varying angles of incidence of the impinging flux. As shown in Brett & Smith (1993), the effect of irradiation on the atmosphere will decrease with increasing angles of incidence. Note that adopting the maximum case, i.e no redistribution at all of $`F_{\mathrm{inc}}`$, affects significantly the outer atmospheric profile, but only slightly the inner profile. We did not consider horizontal energy flows, which may exist in a real planet receiving a varying amount of incident flux over its surface. But our assumption of isotropic incident flux, with the maximum amount of flux allowed, should provide an upper limit for the expected effects of irradiation on the evolution of a planet.

Although present atmosphere models still have some shortcomings, due either to incomplete molecular opacities or to dust treatment (see §2.2), the inner profiles of the irradiated models must be changed and heated drastically to provide the effect required on evolutionary models to match the radius of HD209458b. Fig. 7 indeed shows that the contraction of the irradiated sequence proceeds too rapidly after the first Myr of evolution to maintain a radius close to the observed value. At an age $`t_1`$ = 1 Myr, the model has a radius $`R_1=2R_\mathrm{J}`$ and an intrinsic luminosity $`\mathrm{log}L_1/L_{}=5.12`$, corresponding to $`T_{\mathrm{eff}}`$ = 670K. This corresponds to a characteristic Kelvin-Helmholtz timescale $`\tau _{\mathrm{KH}}=Gm_\mathrm{p}^2/(R_1L_1)10`$Myr. In fact, the model has already reached a radius of $`R_{10}=1.48R_\mathrm{J}`$, about the radius of HD209458b, after only 10 Myr (see Fig. 7). In order to slow down the planet contraction to reach the observed value at $`5`$ Gyr, the inner profiles of the present irradiated atmosphere models, for a given $`T_{\mathrm{eff}}`$, should be drastically modified. Such a major modification seems unlikely, for a fixed incident flux. Note that assuming a constant flux from the parent star during the whole planet evolution overestimates the total incident flux received by the planet over $``$ 5 Gyr. For a parent star mass $`1.06M_{}`$ (Cody & Sasselov 2002), most of the star evolution after the first 1 Myr proceeds at a luminosity $`L`$ smaller than its value at 5 Gyr.

Finally, possible uncertainties due to abundance effects, such as non-solar metallicity and/or helium abundance are not expected to affect significantly the present results. To estimate such uncertainties, we have computed a grid of irradiated atmosphere models for an over-solar metallicity \[M/H\] = +0.3 and the corresponding evolutionary sequence for a mass $`m_\mathrm{p}`$ = 0.69 $`M_\mathrm{J}`$. After a few Gyr evolution, the radius of the later sequence is essentially the same as for the solar metallicity case. Given the optimization of irradiation effects provided by our assumptions, and the huge effect required on the inner profile to reduce the mismatch between observed and predicted radii, we do not expect uncertainties in the present models to be the source of the discrepancy.

### 4.2 Observed versus theoretical radius

The definition of the radius in low mass stars, brown dwarfs or isolated giant planets is usually not a matter of confusion, given the negligible extension of their atmosphere compared to the total radius of the object. The picture could be different in the case of irradiated atmospheres, where extension effects due to the large heating of the upper layers may not be negligible (Seager & Sasselov 2000; Hubbard et al. 2001).

As discussed in Baschek et al. (1991), the condition of compactness for a photosphere in hydrostatic equilibrium is $`H_\mathrm{P}/r<<1`$, whit $`H_\mathrm{P}`$ the pressure scale height and $`r`$ the radial distance to the center. This condition is known to be perfectly fulfilled in (non-irradiated) low mass objects (see Chabrier & Baraffe 1997), where the extension of the photosphere is usually less than 1% the total size of the object. Thus the radius can unambiguously be defined, quoting Baschek et al. (1991), as the distance of the atmosphere to the object center. In all our previous work, we fix the boundary condition between atmosphere and inner structure at $`\tau _{\mathrm{std}}`$=100, knowing that $`R(\tau _{\mathrm{std}}=100)`$ is essentially the same as $`R(\tau _{\mathrm{std}}1)`$ (see Chabrier & Baraffe 1997). As already mentioned, we define $`\tau _{\mathrm{std}}`$ at 1.2 $`\mu `$m, which corresponds to the peak of the flux emitted by cool (non-irradiated) objects. Usually, the region where $`\tau _{\mathrm{std}}1`$ is close to the region where $`\tau _{\mathrm{Rosseland}}1`$ (for the present irradiated models as well). Above this region, the atmosphere contains little mass and contributes negligibly to the luminosity. The evolutionary calculations for irradiated models presented in §3.2 determine also the radius at $`\tau _{\mathrm{std}}=100`$. Within the present assumptions of irradiation, the atmospheric extension between $`\tau _{\mathrm{std}}`$=100 and $`\tau _{\mathrm{std}}1`$ represents only 1-2% of the total radius. The presently calculated theoretical radius is thus essentially equivalent to a photospheric radius at 1.2 $`\mu `$m.

However, the observed radius of HD209458b (e.g the one estimated by Cody & Sasselov 2002) is based on the analysis of optical light curves. It corresponds to a region of the atmosphere where the optical depth is near unity at $`0.6\mu `$m, which is near $`\tau _{\mathrm{std}}=10^2`$ in our atmosphere models. Therefore, the radius predicted by the evolutionary models is not equivalent to the observed radius. If the atmospheric extension and the opacity of the atmosphere are large at the observed wavelengths, the measured radius could be different from the radius predicted by the evolutionary calculation. However, based on our irradiated atmosphere models with the gravity predicted by the evolution and the age of HD209458b (i.e. $`\mathrm{log}g`$=3.2), the atmospheric extension between $`\tau _{\mathrm{std}}=100`$ and $`\tau _{\mathrm{std}}=10^2`$ (where the optical depth is close to unity at $`0.6\mu `$ m) is very small, namely 0.05 $`R_\mathrm{J}`$, compared to the overall radius (Barman et al. 2003, in preparation). Adding this value to the radius predicted by the models at 5 Gyr yields an optical-depth radius at $`0.6\mu m`$ of $``$ 1.14 $`R_\mathrm{J}`$, still 22% less than the observed value. For younger planets or planets undergoing stronger irradiation effects, with much lower gravities, the extension will be more important and should be taken into account for a consistent comparison between theoretical and observed radii, as already stressed by Seager & Sasselov (2000) and Hubbard et al. (2001).

In the same vein, Hubbard et al. (2001) estimate a radius of 94 430 km (1.32 $`R_\mathrm{J}`$) at a pressure of 1 bar, based on a detailed analysis of physical effects influencing the observed light curve of HD 209458b. At 5 Gyr, our models predict a radius at 1 bar of $``$ 1.1 $`R_\mathrm{J}`$, 18% smaller than the Hubbard et al. (2001) estimate. Such a discrepancy is consistent with the afore-mentioned mismatch for the radius at $`0.6\mu `$m.

### 4.3 Other sources of energy deposit

If irradiation effects alone do not explain the large observed radius of HD 209458b, other sources of energy must be invoked. Tidal interactions between the star and the planet can provide a source of energy associated to the synchronization and/or circularization of the planet orbit, dissipated within the planet (Lubow et al. 1997; Rieutord & Zahn 1997; Bodenheimer et al. 2001). However, as discussed recently by Guillot & Showman (2002) and Showman & Guillot (2002), these processes are efficient only during the early stages of the planet evolution. Estimates based on the current understanding of such processes yield typical circularisation timescale $`\tau _{\mathrm{circ}}10^8`$ yr (Bodenheimer et al. 2001) and synchronisation timescale $`\tau _{\mathrm{syn}}<10^8`$ yr (Lubow et al. 1997; Rieutord & Zahn 1997). Such an energy source seems thus unlikely to slow down the long term evolution of the planet, unless a second planet orbiting HD209458a is present. Such a detection has been claimed very recently in the literature (Bodenheimer, Laughlin & Lin, 2003) but remains to be confirmed unambiguously. Showman & Guillot (2002) suggested that downward transport of kinetic energy produced by atmospheric circulation could be dissipated in the planet interior, leading to a substantial deposit of energy. Within the present input physics and treatment of irradiation, we can estimate the amount of energy required to reach the radius of HD 209458b. As in Guillot & Showman (2002), we arbitrarily add an extra term of energy generation $`\dot{ϵ}_{\mathrm{extra}}`$ in the energy equation at different depths. We have explored several cases displayed in Fig. 8. We add a total amount $`L_{\mathrm{extra}}=\dot{ϵ}_{\mathrm{extra}}𝑑m`$ in a region of mass $`\mathrm{\Delta }m`$ enclosed between the surface and an arbitrary depth at mass shell $`m_1`$ (i.e $`\mathrm{\Delta }m=m_\mathrm{p}m_1`$). Various tests indicate that an amount of energy $`L_{\mathrm{extra}}10^{27}5\times 10^{27}`$ erg.s-1 dissipated along the internal adiabat yields a radius within the error bars of the observed value (see Fig. 8). As expected, the larger the fraction of $`L_{\mathrm{extra}}`$ deposited in the convective layers, the more important the effect. Note that the case displayed in Fig. 8 with $`L_{\mathrm{extra}}10^{27}`$ erg.s-1 dissipated all over the star ($`\mathrm{\Delta }m/m_\mathrm{p}`$ = 1, dashed curve) is equivalent to depositing the same amount of energy only at the very center. Our quantitative estimates are in general agreement with Guillot & Showman (2002). Such an amount of energy represents more than 100 times the intrinsic luminosity $`L_{int}`$ of the planet, which is $``$ 1025 erg.s-1 at 1 Gyr and $``$ 2 1024 erg.s-1 at 5 Gyr (see Fig. 7). However, it represents only $`1\%`$ of the incoming luminosity, $`L_{inc}=2\pi R_p^2F_{inc}10^{29}`$ erg s-1, which largely dominates the planet energetic balance. Thus, an alternative possibility is the release of an external source of energy caused by the incident radiation (see e.g. Showman & Guillot, 2002). As illustrated in Fig. 8, however, the extra source of energy must be dissipated at the top of the internal adiabat, i.e. at a much deeper level than the penetration of the incident photons ($`>>\tau =1`$).

## 5 Conclusion

We have presented calculations describing the evolution of cool brown dwarfs and extra-solar gaseous planets. The present models reproduce the main trends of observed methane-dwarfs in near-IR color-magnitude diagrams ($`J`$-$`K`$, $`K`$-$`L`$). Problems still remain at wavelengths $`<`$ 1 $`\mu `$m, with a flux excess predicted in the I-bandpass. The treatment of atomic line broadening in such dense objects may be the source of the present discrepancy. The models fail to reproduce the coolest L-dwarfs and a detailed treatment of dust diffusion/sedimentation is required for a correct description of the transition region between L- and T- dwarfs. Work is in progress in this direction.

We have included the effects of irradiation, coupling irradiated atmosphere profile and inner structure, and providing consistent evolutionary models for irradiated planets. The effect of irradiation are shown to modify significantly the mechanical (mass-radius) and thermal evolution of irradiated EGPs. However, a significant discrepancy (26%) remains between the theoretical and observed radii of the transit planet HD209458b. We have explored possible uncertainties inherent in the models to explain such a discrepancy. Although solving these uncertainties may modify the outer structure of the models (extension of the atmosphere, albedo) and perhaps slightly reduce the discrepancy, none of the uncertainties is likely to modify significantly the inner entropy profile of the models, which determines the radius of the planet. Indeed, a drastic modification on the inner thermal structure is required to bring the theoretical radius in agreement with the observed one.

In summary, we do not expect irradiation effects alone to explain the large observed radius of HD209458b. In the same vein, tidal interactions will affect only the early stages of evolution of the planet but will probably be dissipated too rapidly to affect the long term contraction of the object. Other sources of energy, representing about 100 times the intrinsic luminosity of the planet, seem to be required to explain the observed radius. The first extra-solar planet transit thus remains a challenge for theory. Detection of other transits is now crucial to conclude whether HD209458 is a peculiar system, whether a second planetary companion is confirmed or not, or whether we are missing something in the current understanding of close-in giant planets.

Note:  Isochrones for $`t`$ 1 Myr of the COND models (from 0.5 $`M_\mathrm{J}`$ to 0.1 $`M_{}`$) are available at:
http://www.ens-lyon.fr/~ibaraffe/COND03\_models

###### Acknowledgements.

We are very grateful to H. Harris and S. Leggett for providing data under ascii files and to Doug Lin for mentioning the possible detection of a second planet. We thank our anonymous referee for valuable comments. I.B thanks the Max-Planck Institut für Astrophysik in Garching for hospitality during elaboration of part of this work. This research was supported in part by the LTSA grant NAG 5-3435, the NASA EPSCor grant to Wichita State University, NSF grants AST-9720704 and AST-0086246, NASA grants NAG5-8425, NAG5-9222, as well as NASA/JPL grant 961582 to the University of Georgia. This work was supported in part by the Pôle Scientifique de Modélisation Numérique at ENS-Lyon. Some of the calculations presented in this paper were performed on the IBM pSeries 690 of the Norddeutscher Verbund für Hoch- und Höchstleistungsrechnen (HLRN), on the IBM SP “Blue Horizon” of the San Diego Supercomputer Center (SDSC), with support from the National Science Foundation, on the IBM SP and the Cray T3E of the NERSC with support from the DoE, and using the computer facilities at Centre d’Etudes Nucléaires de Grenoble, CINES and IDRIS. We thank all these institutions for a generous allocation of computer time.


text:
# Joint CMB and Weak Lensing Analysis; Physically Motivated Constraints on Cosmological Parameters

## Abstract

We use Cosmic Microwave Background (CMB) observations together with the Red-sequence Cluster Survey (RCS) weak lensing results to derive constraints on a range of cosmological parameters. This particular choice of observations is motivated by their robust physical interpretation and complementarity. Our combined analysis, including a weak nucleosynthesis constraint, yields accurate determinations of a number of parameters including the amplitude of fluctuations $`\sigma _8=0.89\pm 0.05`$ and matter density $`\mathrm{\Omega }_m=0.30\pm 0.03`$. We also find a value for the Hubble parameter of $`H_0=70\pm 3\text{km}\text{s}^1\text{Mpc}^1`$, in good agreement with the Hubble Space Telescope (HST) key–project result. We conclude that the combination of CMB and weak lensing data provides some of the most powerful constraints available in cosmology today.

Cosmology: Cosmic Microwave Background, Weak Lensing, Large Scale Structure, Dark Matter  preprint: CITA-2003-11

The physics behind the anisotropies we see in the microwave background is well studied and understood leshouches . The evolution of the photon distribution function in the tight coupling era and through decoupling is well inside the linear perturbation regime and is the reason for the CMB’s unique status as a probe of cosmological models. The physical interpretation of the angular power spectrum of primary CMB anisotropies is unambiguous when restricted to the inflationary paradigm and given a suitably parametrized spectrum of initial perturbations.

The recently released WMAP first year results Bennett:2003bz  have revealed the CMB angular power spectrum with unprecedented accuracy to multipoles below $`\mathrm{}=900`$ Hinshaw:2003ex . The results are a stunning confirmation of the acoustic oscillation picture, with perturbations arising from an initial super-horizon spectrum of predominantly adiabatic fluctuations, as predicted for example by simple inflationary models. The measurements of the first two acoustic peaks has confirmed in precise detail earlier detections of the peak/dip pattern on scales below the sound horizon at last scattering toco ; boom ; maxima ; dasi ; cbi ; vsa ; acbar .

On its own, the current picture of the CMB made up of the WMAP results together with high resolution CBI and ACBAR111WMAPext combination observations implies tight constraints on a number of parameters Spergel:2003cb ; Verde:2003ey ; the curvature in units of critical density $`\mathrm{\Omega }_K`$ and various other parameters in the combinations determined by the physical mechanisms which give rise to the observed CMB anisotropy. In addition the measurement of a cross-correlation between the polarization and temperature anisotropy Kogut:2003et ; Page:2003fa  is the first significant detection of reionization in the CMB, which gives a constraint on the optical depth to the last scattering surface.

Although the CMB data alone provide tight constraints on some parameter combinations, other combinations are very poorly constrained due to partial degeneracies. The addition of other data such as measurements of the matter power spectrum $`P(k)`$ is essential to break these degeneracies and tightly constrain the parameters of most interest individually. One way to infer the matter power spectrum is to rely on visible tracers of the (dark) matter distribution such as galaxy redshift surveys or observations of the Lyman-$`\alpha `$ forest. The Lyman–$`\alpha `$ forest gives a way to measure the linear power spectrum of neutral gas at redshifts higher than those probed by galaxy surveys.

The combination of CMB, 2dFGRS Colless:2001gk , and Lyman–$`\alpha `$ forest data croft2002  yields tight constraints on the density of dark matter and vacuum energy, and also reveal an indication of a running of the scalar spectral index characterized by the parameter $`dn_s/d\mathrm{ln}k`$ Spergel:2003cb . However inferring the matter power spectrum using these techniques involves a heuristic treatment of the relation between the tracers and the dark matter usually referred to as ‘biasing’ kaiser . As we enter the much heralded era of precision observations, such heuristic treatments might limit the accuracy with which parameters can be determined. A direct measurement of the power spectrum would not suffer from such limitations.

In terms of physical interpretation, measurements of the lensing signal induced by the LSS (cosmic shear) hold a unique position in the growing set of observational tools available to cosmologists; it is a direct probe of the projected matter power spectrum over a redshift range determined by the lensed sources and over scales ranging from the linear to non–linear regime. The intervening LSS induces a small, coherent correlation in the shapes of the background galaxies which nowadays can be measured accurately Bacon:2000sy ; Hoekstra:2002cj ; Kaiser:2000if ; Ludo2000 ; vanWaerbeke:2002pw . The use of weak lensing data is not without challenges: the small signal requires large survey areas and a careful removal of the observational distortions. However separation of the shear signal into gradient (“E-Type”) and curl (“B-Type”) components provides a control on systematics including the presence of intrinsic alignments of nearby galaxies or systematically induced distortions in the image. The RCS 53 sq. deg. results used in this work Hoekstra:2002xs  have a low B-Type component on large scales together with a well determined redshift distribution of background sources.

In this letter we present results from cosmological parameter fits using only CMB and weak lensing data. The motivation for this approach is to provide constraints on parameters using only observables with robust physical interpretations.

To evaluate the posterior distribution of the parameters of interest from the data we use an extension of the publicly available Markov Chain Monte Carlo package cosmomc222http://cosmologist.info/cosmomc/, as described in Lewis:2002ah . We calculate the likelihood of each cosmological model with respect to a combination of CMB and RCS data. The CMB data consists of WMAP data below $`\mathrm{}=900`$ and CBI , ACBAR , and VSA band powers above $`\mathrm{}=800`$ where the WMAP data is noise dominated and hence the band powers are essentially independent. To compare each angular power spectrym to the WMAP data we use the likelihood calculation routine made available by the WMAP team333http://lambda.gsfc.nasa.gov/ Verde:2003ey .

For each model we also calculate the mass aperture variance $`M_{ap}^2(\theta )`$ Schneider:1997ge  at each aperture $`\theta `$ sampled by the RCS results Hoekstra:2002xs . The mass aperture variance is a narrow filter of the convergence power spectrum $`P_\kappa (\mathrm{})`$ defined as

$$P_\kappa (\mathrm{})=\frac{9}{4}\left(\frac{H_0}{c}\right)^4\mathrm{\Omega }_m^2_0^{\chi _H}𝑑\chi \frac{g^2(\chi )}{a^2(\chi )}P_{3D}(\frac{\mathrm{}}{f_K(\chi )};\chi ),$$
(1)
where $`\chi `$ is the radial coordinate and $`f_K(\chi )`$ is the comoving angular diameter distance to $`\chi `$. $`P_{3D}(k,\chi (z))`$is the 3D power spectrum of matter fluctuations. For each model we use the matter power spectrum calculated by camb Lewis:1999bs  at $`z=0`$ and rescale to $`z>0`$ using the solution for growth of linear perturbations. To include the non–linear contribution to the power spectrum at each redshift we use the halofit procedureSmith:2002dz . The procedure has been calibrated using numerical simulations of structure formation and is significantly more accurate than the previous procedure by Peacock & Dodds pd96 . In particular it reproduces accurately, with rms errors of a few percent, the full non–linear spectrum in standard $`\mathrm{\Lambda }`$CDM models down to scales $`k10h`$Mpc-1. The accuracy of the halofit procedure is adequate for current weak lensing data although future surveys will require more accurate estimates of the full, non–linear power spectrum. This will most probably require the use of large numbers of numerical simulations to calibrate directly the non–linear evolution in the full parameter space.

The function $`g(\chi )=_\chi ^{\chi _H}𝑑\chi ^{}p(\chi ^{})f_K(\chi ^{}\chi )/f_K(\chi ^{})`$ is the source–averaged distance ratio where $`p(\chi (z))`$ describes the redshift distribution of sources in the shear survey which is approximated by the function $`p(z)(z/z_s)^\alpha \mathrm{exp}\left[(z/z_s)^\beta \right]`$. The values $`\alpha =4.7`$, $`\beta =1.7`$, and $`z_s=0.302`$ give the best fit to the observed redshift distribution. To allow for the uncertainty in the mean redshift of the distribution we marginalize over the range of values $`z_s[0.274,0.337]`$ for each likelihood evaluation. This corresponds to the $`\pm 3\sigma `$ range indicated by the $`\chi ^2`$ of the fit to the photometric redshift distribution. The mean redshift for this choice of parameters is $`z=0.540.66`$. We assume a Gaussian prior for $`z_s`$ in this range.

For each model sampled by the Monte Carlo chain we calculate the log likelihood with respect to the RCS data

$$\mathrm{ln}L=\frac{1}{2}\left(\stackrel{~}{M_{ap}^2}_iM_{ap}^2_i\right)C_{ij}^1\left(\stackrel{~}{M_{ap}^2}_jM_{ap}^2_j\right),$$
(2)
where $`\stackrel{~}{M_{ap}^2}_i`$ is the observed mass aperture variance at an aperture $`\theta _i`$ and $`C_{ij}`$ is the covariance matrix of the data Hoekstra:2002xs . This result is added to the log likelihoods from the CMB fit for the same model to obtain the full likelihood with respect to both CMB and RCS data.

We sample the probabilities with respect to six basic cosmological parameters: the physical densities of baryons $`\mathrm{\Omega }_bh^2`$, and cold dark matter $`\mathrm{\Omega }_ch^2`$, the Hubble parameter $`H_0100h\text{km}\text{s}^1\text{Mpc}^1`$, a reionization redshift parameter $`z_{\text{re}}`$, and a constant spectral index $`n_s`$ and amplitude $`A_s`$ of the initial scalar curvature perturbations. We assume the universe is spatially flat, with purely adiabatic perturbations evolving according to General Relativity. The density of a cosmological constant type component $`\mathrm{\Omega }_\mathrm{\Lambda }`$ follows from $`\mathrm{\Omega }_\mathrm{\Lambda }=1.0\mathrm{\Omega }_m`$. We generated sixteen converged Monte Carlo chains using the CMB data only, removed burn in and thinned to obtain fairly independent samples. The matter power spectrum and RCS likelihood was then computed for each sample, and importance sampling used to adjust the chain weights accordingly (see Lewis:2002ah ). The resultant set of weighted samples for the full posterior distribution from the CMB and RCS data were then used to compute our results. The only external prior assumed is a conservative big bang nucleosynthesis (BBN) Gaussian prior of the form $`\mathrm{\Omega }_bh^2=0.022\pm 0.002`$ (1$`\sigma `$) Tytler:2000qf . We include this prior to partially break the remaining $`n_s`$$`\mathrm{\Omega }_bh^2`$$`\tau `$$`A_s`$ degeneracy in the CMB data. The action of this is similar to the $`\tau <0.3`$ prior adopted in the WMAP analysis Spergel:2003cb ; Verde:2003ey .

From the set of samples it is simple to also compute the posterior distribution of other derivable quantities such as the rms amplitude of matter fluctuations on $`8h^1`$Mpc scales assuming linear evolution, $`\sigma _8`$, the total matter density, $`\mathrm{\Omega }_m`$, the optical depth to last scattering, $`\tau `$, and the age of the universe. In this letter we do not consider tensor perturbations, dynamical dark energy candidates, or a running spectral index. We will explore the constraints on these generalized models from CMB and weak lensing data future work.

The set of samples from the full six dimensional parameter space can be used to evaluate marginalized parameter distributions by evaluating the weighted number density of samples ignoring the values of the parameters marginalized over. In Fig. 1 we show the one dimensional marginalized distributions for a number of parameters. Each panel compares the distribution obtained using CMB data with that obtained using CMB and RCS data together; both also include the weak BBN prior discussed above. The effect of adding the weak lensing results is clearly seen in a number of parameters.

In Table 1 we summarize the marginalized constraints for a number of fundamental and derived parameters. We show the results obtained with and without inclusion of the RCS data. The addition of RCS data reduces the errors on $`\sigma _8`$, $`\mathrm{\Omega }_m`$, $`H_0`$, $`\mathrm{\Omega }_\mathrm{\Lambda }`$, and $`\mathrm{\Omega }_ch^2`$. We also show constraints on the ‘classical’ combinations probed by LSS data, namely, the constrained direction $`\sigma _8\mathrm{\Omega }_m^{0.5}`$ and the shape parameter $`\mathrm{\Gamma }\mathrm{\Omega }_mh`$.

It is instructive to look at the marginalized, two dimensional likelihood in the ($`\mathrm{\Omega }_m,\sigma _8`$) plane to understand how drastic improvements in the determination of the two parameters are obtained (Fig. 2). The RCS data alone is near degenerate in a particular direction while CMB data alone provides broad constraints in a quasi–orthogonal direction to RCS. The combination of the two data sets give a much tighter confidence region. The region of intersection in six dimensions has slightly above average CMB likelihood, as is readily assessed using the importance weighted samples, so the data sets are highly consistent even in the full parameter space. The spread of the CMB posterior in the direction of the RCS degeneracy is largely due to the uncertainty remaining in the optical depth, as is clear for the tight constraint for $`\sigma _8e^\tau `$ given in Table 1.

Overall our results are consistent with similar constraints from a combination of CMB, 2dFGRS, and Lyman–$`\alpha `$ data Bennett:2003bz ; Spergel:2003cb  with similar or smaller errors. The values obtained for $`\sigma _8`$ using the WMAP data are somewhat higher than those obtained previously from CMB data due to the new evidence for a significant optical depth and a slightly higher anisotropy amplitude than previous observations indicated Hinshaw:2003ex . This is still lower, although not inconsistent, with estimates of $`\sigma _8`$ from a possible Sunyaev-Zeldovich Effect (SZE) contribution to the CMB power spectrum at high-$`\mathrm{}`$. The latest estimates using the CBI deep–field results Mason:2002tm ; Bond:2002tp  and ACBAR and BIMA Dawson:2002dg  data suggests a value of $`\sigma _8^{SZ}=0.98_{0.21}^{+0.12}`$ Goldstein:2002gf  with large errors due mainly to the non-Gaussian nature of the SZE. Increasingly accurate measurements of the CMB power spectrum at high-$`\mathrm{}`$ will reduce these errors drastically and comparing the two independent determinations of $`\sigma _8`$ will be useful in increasing our understanding of the cluster properties that determine the SZE.

Our result for the Hubble parameter is consistent within 1$`\sigma `$ with the HST key–project result Freedman:2000cf  but has smaller errors. Similarly, the value for the matter density $`\mathrm{\Omega }_m=0.30\pm 0.03`$ is consistent with other determinations. The addition of RCS data leave estimates of the scalar spectral tilt $`n_s`$ essentially unaffected. This is due to the small range of scales probed by the RCS weak lensing results. Future surveys will most certainly have much more leverage on $`n_s`$ as they will probe a range in scales an order of magnitude larger.

We have shown how CMB and weak lensing results can be combined to obtain constraints on cosmological parameters that depend on observations that have simple physical interpretations. Although only first generation weak lensing data are currently available our approach yields results with errors comparable to or even smaller than those obtained using CMB in combination with other types of surveys. These results are encouraging for the use of next generation weak lensing surveys in deriving robust parameter fits. In particular the Canada–France–Hawaii–Telescope (CFHT) Legacy Survey $`170`$ sq. deg. cosmic shear project will be a major step forward in the field of weak lensing.

The increasing accuracy in the determination of the source redshift distribution in future surveys will also help in reducing uncertainties and systematics tied to any intrinsic alignment in the ellipticity of nearby sources. It will also introduce the possibility of resolving separate redshift contributions to the convergence power spectrum (Eq. 1) thus enhancing the parameter fitting ability of the observations.

We conclude that the combined CMB, weak lensing approach to parameter determination already constitutes a competitive alternative to other combinations and holds much promise for future investigations.

It is a pleasure to thank Dick Bond, Ue-Li Pen, and Dmitry Pogosyan for useful discussions. We acknowledge the RCS team for use of their data. Research at CITA is supported by NSERC and the Canadian Institute for Advanced Research. The computational facilities at CITA are funded by the Canadian Fund for Innovation. CRC acknowledges the Kavli Institute for Theoretical Physics (National Science Foundation Grant No. PHY99-07949) where part of this work was carried out.