text:
# Extremely Close-In Giant Planets from Tidal Capture

## 1. Introduction

The discovery of an extrasolar planetary companion orbiting a mere $`0.05\mathrm{AU}`$ from its host star 51 Peg (Mayor & Queloz, 1995) originally came as a surprise, given the seeming impossibility of forming a giant planet so near to its parent star. However, the theorists quickly recovered, and, drawing on older works which studied interactions between protoplanets and their natal disks (Goldreich & Tremaine, 1980; Ward, 1986), invoked planetary migration as the mechanism for delivering the companion to 51 Peg from its birthplace to its current observed position (e.g., Lin, Bodenheimer, & Richardson 1996).

Since this time, considerable progress in the theory of planetary migration, combined with the detection of a large number of additional giant planets on short-period orbits, has left little doubt that this process is important in shaping the architecture of planetary systems. However, many questions remain. For example, the mechanisms by which planetary migration is stopped are uncertain (Lin, Bodenheimer, & Richardson, 1996; Trilling et al., 1998; Kuchner & Lecar, 2002). These mechanisms have the additional burden that they should explain both the observed ‘pile-up’ of planetary companions at periods of $`P=3\mathrm{d}`$ (i.e. Kuchner & Lecar 2002), and the newest discovery of a planet at $`P=1.2\mathrm{d}`$ (Konacki et al., 2003; Sasselov, 2003). In addition, it is clear from examples of extrasolar giant planets orbiting close to their supposed birthplaces, as well as our own solar system, that orbital migration is not always so efficient. A satisfactory explanation for the causes of these variations in the efficiency of migration is lacking.

It seems plausible that many, if not the majority, of disk stars were initially formed in stellar systems – loose associations or open clusters that have since dissolved. Evidence for this comes in part from studies of the nearby moving groups (de Zeeuw et al., 1999). Planets orbiting stars in stellar systems are subject to numerous effects resulting from the dense stellar environment. Planets orbiting at larger distances from their parent stars are generally more sensitive to these effects. Therefore, the structure of planetary systems will be the result of a complex interplay between the local properties of the system that drive planetary formation and migration, and the non-local effects that arise from the star’s environment. Thus, if it is indeed the case that most stars formed in stellar clusters, planetary systems cannot be understood as isolated systems: a complete picture of planetary formation and evolution must consider both local and non-local effects.

Several authors have considered the effects of dense stellar environments on the formation and survival of planetary systems (Armitage, 2000; Davies & Sigurdsson, 2001; Smith & Bonnell, 2001; Bonnell, Smith, Davies, & Horne, 2001; Hurley & Shara, 2002). From these studies it is clear that, if distant planets can form in stellar clusters, they are likely to be stripped from their parent stars, resulting in a free-floating population of planets that will slowly evaporate from the system. Recently, Bonnell et al. (2003) suggested that freely-floating brown dwarfs (BD) in globular clusters will occasionally be tidally captured, resulting in a detectable population of extremely close-in BD companions to main-sequence stars if the frequency of freely-floating BDs is large. Similar considerations can also be applied to planets in stellar systems.

Here I consider the formation and stability of extremely close-in giant planetary companions to stars via tidal capture of free-floating planets. In §2, I summarize the dynamical processes of ionization, evaporation, capture, and tidal decay, and their relevant timescales. In §3, I give a crude estimate for the frequency of extremely close-in giant planets as a function of the parameters of the stellar system. In §4, I consider tidally-captured BD and planetary companions to stars in globular clusters, and interpret the null result of the 47 Tuc transit search (Gilliland et al., 2000). Finally, in §5, I discuss possible implications and prospects for the detection of such a population of planetary companions.

## 2. Dynamical Processes

It is not clear if planets can form around stars in stellar systems at all. In particular, in rich clusters with $`N10^3`$ members, high-mass stars can generate an ultraviolet radiation field that is sufficiently intense to photoevaporate protoplanetary discs in a few hundred thousand years (Armitage, 2000). However, planet formation might still be possible if there exists a substantial delay between high and low-mass star formation. This fact, combined with our incomplete understanding of the efficiency of planetary migration, make it difficult to predict a priori the frequency and distribution of planets orbiting stars in stellar systems. For the purposes of discussion, I will sweep all of these uncertainties under the rug, and simply consider the evolution and consequences of a substantial population of long-period planets around stars in stellar systems.

Binaries in stellar systems evolve under random encounters with other stars in the system. This evolution proceeds according to Heggie’s law (Heggie, 1975), such that soft binaries tend to get disrupted by encounters. A binary is soft when its binding energy is equal to the mean kinetic energy of the stars in the system. Thus there is a critical separation above which binaries will tend to be disrupted,

$$a_{\mathrm{dis}}=\frac{Gm_pM_{}}{2M\sigma ^2},$$
(1)
where $`m_p`$ is the mass of the secondary, $`M_{}`$ is the mass of the primary, $`M`$ is the average mass and $`\sigma `$ is the velocity dispersion of the cluster stars, respectively. Planetary companions to solar-type stars ($`m_pM_{\mathrm{Jup}}`$ and $`M_{}M_{}`$) in typical open clusters ($`M0.5M_{}`$ and $`\sigma 1.5\mathrm{km}\mathrm{s}^1`$ are will tend to get disrupted if $`a0.4\mathrm{AU}`$, or $`P100\mathrm{d}`$.

The timescale for disruption of a planet is approximately the time it takes for the random encounters to change the planet’s energy by an amount equal to its binding energy (Binney & Tremaine, 1987),

$$t_{\mathrm{dis}}=(1+q)\frac{M_{}}{M}\frac{\sigma }{16\sqrt{\pi }GM\nu a\mathrm{ln}\mathrm{\Lambda }}.$$
(2)
Here $`qm_p/M_{}`$, $`\nu `$ is the stellar number density, and $`\mathrm{ln}\mathrm{\Lambda }\mathrm{ln}0.4N`$ is the Coulomb logarithm. The number density and velocity dispersion within a stellar system can vary by several orders of magnitude. For simplicity, here I will adopt these properties averaged over the region interior to $`r_{hm}`$, the half-mass radius of the cluster. The mean stellar density in this area is simply

$$\nu =\frac{\frac{1}{2}NM}{\frac{4}{3}\pi r_{hm}^3}.$$
(3)
For a wide range of equilibrium stellar systems, the half-mass radius can be related to the total mass and velocity dispersion of the system by (Spitzer, 1987),

$$r_{hm}=\frac{0.4GMN}{3\sigma }.$$
(4)
I will adopt this relation throughout unless stated otherwise. For an average stellar density $`\nu 10^3\mathrm{pc}^3`$ and velocity dispersion $`\sigma 1.5\mathrm{km}\mathrm{s}^1`$, the above relation yields $`N7600`$. For the stellar parameters given above, $`t_{\mathrm{dis}}0.5\mathrm{Gyr}(a/2.5\mathrm{AU})^1`$. Thus most planets with $`a2.5\mathrm{AU}`$ in an open cluster will be disrupted over the cluster’s lifetime ($`0.5\mathrm{Gyr}`$).

Since planets are typically liberated due to the accumulated effects of distant encounters, free-floating planets will initially have velocities similar to the stars, and will not escape the system immediately (Hurley & Shara, 2002). However, equipartition of energy will slowly drive the velocity dispersion of free-floating planets to the point where many of the planets have sufficient velocities to evaporate from the system entirely. Assuming that the stellar component of the system dominates the mass, the characteristic timescale for this process is the relaxation timescale (Spitzer, 1987; Binney & Tremaine, 1987),

$$t_{\mathrm{rlx}}=0.34\frac{\sigma ^3}{\nu (GM)^2\mathrm{ln}\mathrm{\Lambda }}.$$
(5)
For the parameters adopted above, this yields $`t_{\mathrm{rlx}}0.03\mathrm{Gyr}`$.

Freely-floating planets will occasionally pass sufficiently close to another star to raise a significant tide. If the energy required to raise this tide is larger than the relative energy of the planet and star at infinity, then the encounter will lead to a bound system (Fabian, Pringle, & Rees, 1975; Press & Teukolsky, 1977; Lee & Ostriker, 1986). This leads to the requirement that the planet must pass within a minimum distance $`a_{\mathrm{cap}}`$ to be captured (Fabian, Pringle, & Rees, 1975),

$$a_{\mathrm{cap}}R_{}\left[\frac{GM_{}}{R_{}\sigma ^2}q(1+q)\right]^{1/6}.$$
(6)
For $`q=10^3`$, $`M_{}M_{}`$, $`R_{}R_{}`$, $`a_{\mathrm{cap}}/R_{}2`$. Thus the planet must pass within two stellar radii to be captured.

The timescale for tidal capture is,

$$t_{\mathrm{enc}}\frac{\sigma }{8\sqrt{\pi }GM_{}\nu a_{\mathrm{cap}}},$$
(7)
where I have assumed that $`GM_{}/2\sigma ^2a_{\mathrm{dis}}/qa_{\mathrm{cap}}`$, which is valid for the cases considered here. I find $`t_{\mathrm{enc}}500\mathrm{Gyr}(a_{\mathrm{cap}}/2R_{})^1`$. Assuming every star forms $`N_p`$ planets that are liberated from their parent stars essentially immediately, the fraction of stars with tidally captured planet can be crudely estimated as $`N_p(t_{\mathrm{rlx}}/t_{\mathrm{enc}})0.01\%N_p`$.

Planets passing within $`a_{\mathrm{coll}}=R_p+R_{}`$ of a star will physically collide with it, resulting in a merger. This therefore sets an absolute lower limit on the period of a tidally captured planet. For gaseous objects supported only by degeneracy pressure (including giant planets and brown dwarfs), $`R_p0.1R_{}`$, roughly independent of mass. For $`R_{}R_{}`$, $`a_{\mathrm{coll}}0.005\mathrm{AU}`$, or a period of $`P0.13\mathrm{d}`$. This also sets a lower limit on the mass of a planet that can be captured, as very low-mas planets must pass within the collision separation in order to raise a sufficient tide to be captured. This limit is found by equating $`a_{\mathrm{coll}}`$ with $`a_{\mathrm{cap}}`$:

$$q\frac{R_{}\sigma ^2}{GM_{}}\left(1+\frac{R_p}{R_{}}\right)^6,$$
(8)
where I have assumed $`q1`$. For the fiducial parameters, this is $`q2\times 10^5`$.

Low-mass planets may also be tidally captured sufficiently close to their parent stars that their radii exceed the Roche limit. The limiting separation for Roche lobe overflow is,

$$a_\mathrm{R}=\left(\frac{3}{q}\right)^{1/3}R_p1.44\left(\frac{q}{10^3}\right)^{1/3}R_{},$$
(9)
where, for the last relation, I have assumed $`R_p=0.1R_{}`$. The Roche limit separation exceeds $`a_{\mathrm{coll}}`$ for $`q2\times 10^3`$. Planets that are captured within their Roche limit will lose mass to their parent star. This mass loss is accompanied by a transfer of angular momentum, resulting in an outward migration of the planet. Thus a planet initially captured inside its Roche limit can still result in a stable system. The lower the mass of the planet, the farther it must migrate to halt mass transfer. Very low mass captured planets may not survive at all.

Planets captured into very close orbits will quickly circularize, on a timescale of $`10^5`$ years (Rasio et al., 1996), and thereafter be subject to tidal decay. The timescale for tidal decay is (e.g. Rasio et al. 1996),

$$t_{\mathrm{dec}}=\frac{f}{t_\mathrm{c}}\frac{M_{\mathrm{cz}}}{M_{}}q(1+q)\left(\frac{R_{}}{a}\right)^8,$$
(10)
where $`M_{\mathrm{cz}}`$ is the mass of the convective zone; for solar-type stars, $`M_{\mathrm{cz}}0.02M_{}`$. Here $`t_\mathrm{c}`$ is the convective eddy turnover timescale, which I will approximate as $`t_\mathrm{c}[M_{\mathrm{cz}}0.2R_{}^2/3L_{}]^{1/3}`$, where $`L_{}`$ is the luminosity of the star (see, e.g., Rasio et al. 1996). For solar-type stars, $`t_\mathrm{c}20\mathrm{d}`$. Thus for tidally-captured planets, $`Pt_\mathrm{c}`$, and the largest eddies can no longer contribute to the total viscosity, resulting is a suppression of the viscous dissipation. Therefore $`f`$ is less than unity. However, the correct form for $`f`$ remains controversial (see Goodman & Oh (1997) for a discussion). Following Sasselov (2003), I will consider both linear suppression of the eddy viscosity (Zahn, 1989),

$$f_L=\left(\frac{P}{2t_\mathrm{c}}\right),Pt_\mathrm{c},$$
(11)
as well as quadratic suppression (Goldreich & Keeley, 1977),

$$f_Q=\left(\frac{P}{2\pi t_\mathrm{c}}\right)^2,Pt_\mathrm{c}.$$
(12)
Although quadratic suppression is theoretically better motivated, linear suppression appears to be in better agreement with the observed timescale of tidal circularization in close binaries (Goodman & Oh, 1997).

For tidally-captured planets with $`a_{\mathrm{cap}}2R_{}`$, the decay timescale for linear suppression is $`t_{\mathrm{dec}}0.1\mathrm{Gyr}`$, whereas for quadratic suppression, $`t_{\mathrm{dec}}80\mathrm{Gyr}`$. Therefore the choice of prescription has an enormous effect on the number of surviving tidally-captured planets: one expects a negligible number of surviving planets for linear suppression, whereas planets are stable over much longer than $`\mathrm{Gyr}`$ for quadratic suppression.

The angular momentum lost from the tidal decay of the planet’s orbit will be transfered to the star, spinning it up. Assuming the star is initially slowly rotating, spin-up of the star will occur on a timescale (Rasio et al., 1996),

$$t_{\mathrm{su}}\frac{I_{}}{I_{orb}}t_{\mathrm{dec}}=\frac{k^2}{q}\left(\frac{R_{}}{a}\right)^2t_{\mathrm{dec}},$$
(13)
where $`I_{}`$ and $`I_{orb}`$ are the moments of inertia of the star and orbit, and $`k^2I_{}/M_{}R_{}^2`$. For $`a=a_{\mathrm{cap}}2R_{}`$, and $`k^2=0.08`$, $`t_{\mathrm{su}}20(10^3/q)t_{\mathrm{dec}}`$. Thus companions with $`q0.02`$ will decay before they spin-up the star sufficiently to synchronize the star’s spin period with the orbital period. Conversely, companions with $`q0.02`$ will synchronize before they decay completely. This equilibrium is stable only if the orbital angular momentum of the planet is more than three times the spin angular momentum of the star (Hut, 1980). This leads to a minimum mass ratio for stability,

$$q3k^2\left(\frac{R_{}}{a}\right)^2.$$
(14)
Thus companions with $`q0.06`$ cannot reach a stable equilibrium, and will always decay. More massive companions may reach a stable equilibrium.

The question of orbital stability via synchronization is further complicated by the fact that the star will likely be spun-up as a result of the capture process itself111I am indebted to I. Bonnell for pointing this out., and thus may be fairly rapidly rotating. This will not only affect the viscous dissipation mechanisms, but will also generally shorten the timescale for synchronization. Given the uncertainties involved in determining the effects of these various process, and since I am primarily interested in planetary-mass companions, for which synchronization is likely less important, I will furthermore ignore synchronization entirely, and assume that all orbits decay on the timescale given by eq. (10). However, I note that this issue has direct bearing on the discussion of the stability of tidally-captured BD companions to stars in 47 Tuc in §4.

Figure 1 summarizes the regions of parameter space in the companion mass-period plane where ionization, tidal capture, and tidal decay are important, as discussed in this section. For this figure, I have assumed as before $`M_{}=M_{}`$, $`M=0.5M_{}`$, $`R_{}=R_{}`$, $`\sigma =1.5\mathrm{km}\mathrm{s}^1`$, and $`\nu =10^3\mathrm{pc}^3`$. In addition, I have assumed a fixed companion radius of $`R_p=R_{\mathrm{Jup}}`$, and a system lifetime of $`0.5\mathrm{Gyr}`$. Assuming quadratic suppression prescription of tidal decay, there is a small sliver of parameter space in the period-mass plane of where planets can be tidally captured and are stable to tidal decay over typical open cluster lifetimes. For linear suppression, there is no region of parameter space where tidally-captured planets are expected to be stable on such timescales. Also shown are known extrasolar planets. Clearly tidally-captured planets occupy a region of parameter space that is relatively disjoint from known companions, and would likely to easily identified as the result of capture (rather than migration).

## 3. Frequency Estimates

The frequency of tidally-captured planets in a stellar system of a given age depends on the competing effects of ionization, evaporation, capture, and orbital decay. The exact number will depend on the balance between the timescales of these various dynamical effects. This will in turn depend on the dynamical properties of the cluster at all times in its evolution, and thus the processes of relaxation, mass segregation, mass loss, etc., should all be considered. Ideally, the most robust way to accomplish this is through detailed N-body simulations. However, such a study is outside of the scope of this paper. Furthermore, such a detailed study is perhaps not warranted, as the number of tidally captured planets will also depend critically on the initial frequency and distribution of planets formed around stars in the cluster. Given that a general theory of planet formation is lacking, this can only be a wild guess at best. I therefore simply provide only a crude estimate for the frequency of tidally-captured planets. This estimate should be good to an order-of-magnitude, and should serve to elucidate the dependence of the frequency on the gross parameters of the stellar system and the input assumptions.

I consider only average properties of the stars and planets in the stellar system. The dynamical properties of the stellar system are specified by the average number density $`\nu `$ and velocity dispersion $`\sigma `$ inside the half-mass radius, and the average mass of the stars, $`M`$. The total number is then set by equations (3) and (4). Unless otherwise stated, I assume $`M=0.5M_{}`$, which is roughly appropriate for mass functions observed in the field, young clusters, open clusters, and globular clusters (see Chabrier 2003). Here and throughout, I adopt radii and bolometric luminosities from Allen (1976) and convective zone masses from Pinsonneault, DePoy, & Coffee (2001).

For my fiducial calculations, I assume that every star in the system originally has four companions with mass $`M_p=M_{\mathrm{Jup}}`$, with a distribution that is uniform in $`\mathrm{log}a`$ between $`550\mathrm{A}\mathrm{U}`$. This choice is primarily motivated by the distribution of massive planets in our own solar system, but is not inconsistent with the distribution of extrasolar planets detected via radial velocity surveys. These surveys are only just becoming sensitive to planets around solar-type stars at $`5\mathrm{A}\mathrm{U}`$; however extrapolations based on current samples indicate that the fraction of stars with Jupiter-mass planets at periods longer than this may be quite large (e.g., Lineweaver & Grether 2003). For the same distribution of $`a`$ but other choices for the frequency of massive planets, one can simply scale all results by $`(N_p/4)`$, where $`N_p`$ is the average number of planets per star.

Free-floating planets in stellar stellar are depleted via evaporation and replenished by disruption of bound planets222Free-floating planets are also depleted via captures, but the rate of captures is generally so small that this can safely be ignored.. The rate of change in the number $`N_{\mathrm{ff}}`$ of free-floating planets per star is thus

$$\frac{\mathrm{d}N_{\mathrm{ff}}}{\mathrm{d}t}=\frac{\mathrm{d}N_{\mathrm{dis}}}{\mathrm{d}t}\frac{\mathrm{d}N_{\mathrm{evap}}}{\mathrm{d}t}.$$
(15)
Here $`\mathrm{d}N_{\mathrm{dis}}/\mathrm{d}t`$ is the rate at which planets are ionized,

$$\frac{\mathrm{d}N_{\mathrm{dis}}}{\mathrm{d}t}=\frac{\mathrm{d}}{\mathrm{d}t}\left[_0^{\mathrm{}}\frac{\mathrm{d}N_p}{\mathrm{d}\mathrm{log}a}𝒫_{\mathrm{dis}}(a)\mathrm{d}\mathrm{log}a\right].$$
(16)
Here $`\mathrm{d}N_p/\mathrm{d}\mathrm{log}a`$ is the initial distribution of planetary companions, which is constant between $`550\mathrm{AU}`$ and zero otherwise, and $`𝒫_{\mathrm{dis}}1e^{t/t_{\mathrm{dis}}(a)}`$ is the probability that planet with separation $`a`$ will be disrupted. The rate at which planets are evaporated is

$$\frac{\mathrm{d}N_{\mathrm{evap}}}{\mathrm{d}t}=\frac{\xi _eN_{\mathrm{ff}}}{t_{\mathrm{rlx}}},$$
(17)
where $`\xi _e=0.156`$ is the evaporation probability appropriate for test particles (Spitzer, 1987). The number of free-floating planets as a function of time can be found by integrating eq. (15) from $`t=0`$.

Tidally-captured planets are depleted via orbital decay and replenished via new captures of free-floating planets. The rate of change in the number $`N_{\mathrm{sur}}`$ of tidally-captured planets per star is thus

$$\frac{\mathrm{d}N_{\mathrm{sur}}}{\mathrm{d}t}=\frac{\mathrm{d}N_{\mathrm{cap}}}{\mathrm{d}t}\frac{\mathrm{d}N_{\mathrm{dec}}}{\mathrm{d}t}=\frac{N_{\mathrm{ff}}}{t_{\mathrm{enc}}}\frac{N_{\mathrm{sur}}}{t_{\mathrm{dec}}}.$$
(18)
I determine the mean capture cross section by averaging the width of the region between $`a_{\mathrm{cap}}`$ and $`a_{\mathrm{coll}}`$, which depends on the relative velocity between the planet and star, over a Maxwellian (relative) velocity distribution with dispersion $`\sigma `$. I then evaluate $`t_{\mathrm{enc}}`$ for this cross section, and $`t_{\mathrm{dec}}`$ for the average semi-major axis of the captured planets. The number of surviving free-floating planets is then found by integrating eq. (18).

Note that I have ignored the effects of mass segregation. Mass segregation will decrease the rate at which planets are captured, as it will partially decouple the dynamical interactions between the stellar and free-floating planetary populations. This will generally be a relatively small correction, as most planets are tidally captured during the first few relaxation times of the cluster.

Figure 2 shows the resulting evolution of the planetary population for three types of stellar systems, characterized by the parameters $`\nu `$ and $`\sigma `$. These are $`\nu =10^4\mathrm{pc}^3`$ and $`\sigma =10\mathrm{km}\mathrm{s}^1`$, appropriate to a globular cluster, $`\nu =10^3\mathrm{pc}^3`$, $`\sigma =1.5\mathrm{km}\mathrm{s}^1`$, appropriate to a rich open cluster, and $`\nu =10^2\mathrm{pc}^3`$, $`\sigma =0.6\mathrm{km}\mathrm{s}^1`$, appropriate to a loose association. These systems have a total number of stars and half-mass radii equal to $`N7\times 10^5`$, $`r_{hm}=2\mathrm{pc}`$ (globular cluster), $`N7600`$, $`r_{hm}1\mathrm{pc}`$ (rich open cluster) and $`N1500`$, $`r_{hm}1.2\mathrm{pc}`$ (association). These values are not meant to be definitive, but merely span the interesting range of possible systems. Note that $`r_{hm}`$ does not enter into the calculations directly, and $`N`$ enters only logarithmically via the Coulomb logarithm.

Due to their high stellar density and velocity dispersions, globular clusters generally have relatively small disruption timescales and long relaxation timescales. Therefore, planets are liberated from their parent stars very quickly, and evaporation is slow, making the dynamics particularly simple: the number density of free-floating planets is roughly constant, and $`\mathrm{d}\mathrm{log}N_{\mathrm{sur}}/\mathrm{d}\mathrm{log}t`$ is approximately constant when tidal decay is negligible. At the typical globular cluster lifetime of $`10\mathrm{Gyr}`$, $`N_{\mathrm{sur}}0.1\%`$.

The evolution of rich open clusters is more complicated due to their short relaxation times. From inspection of Figure 2, it is clear that the majority of planets are tidally captured within a time $`t_{\mathrm{rlx}}/\xi _e0.2\mathrm{Gyr}`$, after which the number of surviving planets is approximately constant at $`N_{\mathrm{sur}}0.03\%`$. The lifetime of open clusters is uncertain, but rich clusters have been observed with ages of a $`\mathrm{Gyr}`$ or more (Kalirai et al., 2001; Burke et al., 2003a). Less massive clusters are probably dissolved on shorter timescales.

Loose associations have even shorter relaxation times, but are otherwise similar to open clusters. However one important difference is that, because of their smaller velocity dispersions, planets are captured at larger separations. Due to the extremely strong dependence of the tidal decay rate on the semi-major axis, this has a profound effect on the survival of tidally-captured planets for the linear prescription for the suppression of eddy viscosity: even under such a prescription, one expects a significant number of surviving tidally-captured planets at a $`\mathrm{Gyr}`$ or less. Loose associations are likely to disperse fairly rapidly, on timescales of a few hundred million years or less. At $`10^8`$ years, $`N_{\mathrm{sur}}0.005\%`$

## 4. Brown Dwarfs and Planets in 47 Tuc

Because of their high encounter rates and large numbers of stars, globular clusters are expected to contain a large number of systems formed via tidal capture. Using the results from the previous section, a system with $`\nu 10^4\mathrm{pc}^3`$ and $`\sigma 10\mathrm{km}\mathrm{s}^1`$ should have $`NN_{\mathrm{sur}}=7\times 10^5\times 0.001700`$ systems of tidally-captured planets. This process is not limited to planets stripped from their parent stars, of course: any planets and brown dwarfs formed in isolation will also be tidally captured by stars in globular clusters. The initial mass function of single objects is observed to be approximately flat down to well below the hydrogen burning limit in many young clusters (Chabrier 2003, but see Briceño et al. 2002). Such a mass function implies a similar number density of stars and brown dwarfs. If this mass function is universal then globular clusters should also contain a large number of stars with tidally-captured BD companions (Bonnell et al., 2003).

Because of their extremely close orbits, tidally-captured planetary and BD companions have a very high probability ($`25\%`$) of transiting their parent star, making photometric monitoring an ideal way of detecting such companions. Gilliland et al. (2000) monitored the globular cluster 47 Tuc continuously for $`8\mathrm{d}`$ using the Hubble Space Telescope, but found no transits. This null result can be used to place constraints on the initial frequency of planetary companions as well as freely-floating planets and BDs formed in isolation (Bonnell et al., 2003).

Figure 3 shows the expected number of tidally-captured companions per star in 47 Tuc as a function of the mass of the companion, for both linear and quadratic prescriptions for viscosity suppression. I have adopted parameters appropriate to 47 Tuc, namely $`\nu =10^4\mathrm{pc}^3`$, $`\sigma =12\mathrm{km}\mathrm{s}^1`$, $`r_{hm}=3.9\mathrm{pc}`$, $`M=0.56M_{}`$ and an age of $`13\mathrm{Gyr}`$ (Djorgovski, 1993; Gebhardt, Pryor, Williams, & Hesser, 1995; Paresce & De Marchi, 2000). Note that here I do not assume eq. (4). Figure 3 assumes an initial frequency of one object (planet or brown dwarf) per star. Due to the extremely high disruption rate, these results are essentially independent of whether the objects are initially free-floating, or bound on orbits with $`a1\mathrm{AU}`$. For the quadratic prescription for viscosity suppression, the expected frequency of tidally-captured BD ($`0.05M_{}`$) companions is $`N_{\mathrm{sur}}0.002`$ companions per star, whereas the frequency of tidally-captured Jupiter-mass planets is $`N_{\mathrm{sur}}0.0002`$.

Gilliland et al. (2000) claim that, if $`f=0.81.0\%`$ of the 34,091 main-sequence stars they monitored had $`1.3R_{\mathrm{Jup}}`$ companions at periods of $`3.5\mathrm{d}`$, then they should have seen 17 planets. This implies a 95% confidence level (c.l.) upper limit of $`3f/17=0.140.18\%`$ to the fraction of main-sequence stars with such companions. This null result can also combined with the prediction for the number of tidally-captured companions $`N_{\mathrm{sur}}`$ to determine an upper limit to the initial frequency $`N_p`$ of objects per star of

$$N_p\frac{3}{17}\frac{f}{N_{\mathrm{sur}}}\frac{𝒫_T(3.5\mathrm{d})}{𝒫_T(P_{\mathrm{cap}})}.$$
(19)
Here $`𝒫_T=(R_p+R_{})/a`$ is the transit probability. Figure 3 shows the resulting 95% c.l. upper limit to the initial frequency of objects as a function of the mass of the object. The shaded region comes from varying $`f`$ in the range $`0.8\%f1.0\%`$.

The null result implies that, on average, each star in 47 Tuc originally formed $`1`$ planet with mass $`M_pM_{\mathrm{Jup}}`$, or $`0.4`$ planets with $`M_p5M_{\mathrm{Jup}}`$. Less than 50% of stars in 47 Tuc could have formed solar-system analogs. While this limit is considerably weaker than the original limit placed by Gilliland et al. (2000) of $`0.1\%`$, it is important to emphasize that this original limit applies only to planets in close ($`3\mathrm{d}`$) orbits. Since planets in such small orbits were not formed in situ, and must be the result of migration, drawing conclusions about the formation of planets at larger separations based on the lack of planets at small separations is subject to uncertainties in the timescale for orbital migration. If the migration timescale were greater than or of order the disruption timescale, the planets would have been stripped from their parent stars before they had a chance to migrate sufficiently close to be detectable via transits. In contrast, the limit on the frequency of planets implied by the lack of tidally-captured planets, although weaker, is essentially independent of where they were initially formed. In fact, one can conclude the initial number of free-floating planets (planets formed in isolation via gas fragmentation and collapse) is also less than or equal to the number of stars.

Under a similar set of assumptions as adopted here, but using a more sophisticated simulation and neglecting orbital decay, Bonnell et al. (2003) found an upper limit to the BD frequency in 47 Tuc of $`15\%`$. Figure 3 shows the upper limit derived using the formalism in §3, but neglecting orbital decay. I find $`12\%`$, in good agreement with Bonnell et al. (2003). This suggests that, in the absence of orbital decay, the predicted frequency of tidally-captured BD companions should be reasonably robust.

For quadratic suppression of eddy viscosity, the timescale for orbital decay of tidally-captured BD ($`q0.1`$) companions is $`t_{\mathrm{dec}}9\mathrm{Gyr}`$. Therefore, approximately $`\mathrm{exp}(9/13)50\%`$ of the captured BD will still be present today. Thus tidal decay is significant at the current age of 47 Tuc for tidally-captured BD companions. Figure 3 shows the upper limit including orbital decay. The limit is revised upward by $`2`$ relative to the estimate of Bonnell et al. (2003). The lack of transiting BD companions in 47 Tuc implies that the initial frequency of BDs relative to stars is $`25\%`$.

The upper limit of $`25\%`$ on the relative frequency of BDs in 47 Tuc is inconsistent with the observed mass functions of most young clusters and the field, for which BDs and stars exist in equal numbers (Chabrier, 2003), although it is consistent with the frequency of BDs in the Taurus star-forming region (Briceño et al., 2002). At face value, this is an indication that the initial mass function compact objects is not universal (Bonnell et al., 2003). Unfortunately, this conclusion relies heavily on the assumed value for the tidal decay timescale. Different assumptions for the mechanisms for tidal dissipation will give rise to radically different conclusions. For example, in Figure 3 I show the predicted frequency of tidally-captured companions under the assumption of linear suppression of eddy viscosity. The frequencies are always $`3\times 10^6`$, i.e. $`5`$ in the entire cluster. This low number results from the fact that the tidal decay timescale under linear suppression is smaller by a factor of $`\sqrt{f_Q}10^4`$. Thus $`t_{\mathrm{dec}}4\times 10^6\mathrm{yrs}`$, and the frequency of tidally-captured objects is only $`t_{\mathrm{dec}}/t_{\mathrm{enc}}10^6`$. Obviously the resulting upper limits are not very interesting. Even if one accepts that the quadratic suppression of eddy viscosity is correct, there are still some ambiguities. For example, there exist discrepancies in the quoted form for $`f_Q`$: Goodman & Oh (1997) and Sasselov (2003) use the form adopted here, namely $`f_Q=(P/2\pi t_c)^2`$, but Rasio et al. (1996) adopt $`f_Q^{}=f_Q\pi ^2`$, i.e. larger by almost an order of magnitude. This factor translates directly to the implied upper limit on the relative frequency, weakening the constraint considerably (see Fig. 3). Note that this would also weaken the constraint on planets in 47 Tuc, but by a smaller factor than the constraint on BDs.

Thus, given the current uncertainties in the physics of tidal dissipation, I conclude that no robust inferences can made about the frequency of brown dwarf companions in 47 Tuc, and any conclusions regarding the universality of the initial mass function are probably premature.

## 5. Implications and Prospects for Detection

It is clear that, under certain sets of assumptions, extremely close-in tidally-captured massive planets should exist around stars formed in dense stellar systems. What are the prospects for the detection of such companions, and what would be the implications of any such detections?

Given their extremely close orbits, the most promising method for detecting companions is transits. For a planet captured at an orbit of $`a2R_{}`$, the transit probability is $`𝒫_T50\%`$. The duration of the transit is

$$t_T=\frac{P}{\pi }\mathrm{arcsin}\left[\sqrt{𝒫_T^2\mathrm{cos}^2i}\right],$$
(20)
where $`i`$ is the inclination. Thus for a central transit, the duty cycle (fraction of time in transit) is quite large: $`t_T/P15\%`$. For a solar-type primary, the transits will last $`50\mathrm{min}`$ and recur every $`5\mathrm{hours}`$. The transit depth will be $`1\%`$. Given the extremely close proximity of the planet to the star, one might worry that the planet will induce additional photometric variations on the star. The reflected light from the companion will induce an fractional flux variation of order $`\mathrm{\Delta }F/Fp(R_p/a)^22\times 10^3`$, where $`p`$ is the geometric albedo, and I have adopted $`p=2/3`$ as appropriate for Lambert sphere scattering. Although five times lower than the signal from the transit, a signal of this magnitude may be large enough to be detected. In fact, depending on the selection criterion used to choose targets for the transit search in 47 Tuc, this may be cause for concern in interpreting this result. The planet will also induce photometric variations of the star from tidal and Doppler effects, but these will likely be small for planetary companions (Loeb & Gaudi, 2003). They may, however, be significant for more massive companions. Thus, barring an unforeseen effects that induce additional photometric variability, it should be possible to cleanly detect tidally-captured planets via transits.

Because of their short periods and relatively large duty cycles, the observational requirements of a transit search for tidally-captured planets are generally much less severe than transit searches for more distant planets. In particular, complete coverage of several phases can be achieved in less than a few days. This is especially important for ground-based transit searches, as the coherence time of weather patterns is typically several days. This proves to be devastating when searching for planets with periods of order the coherence time, but nearly optimal for very short periods. Furthermore, aliasing is minimized because the period of the planet is typically less than the duration of a typical clear observing night.

Tidally-captured planets will also induce radial velocity variations in the host star, of amplitude $`K340\mathrm{m}\mathrm{s}^1(5\mathrm{hr}/P)(M_p/M_{\mathrm{Jup}})`$. For reasonably bright sources, this is well within reach of current instrumentation. However, the presence of such a close planet may induce additional radial velocity variability, making detection of the planetary signal difficult. Furthermore, the extremely low expected frequency of tidally-captured planets makes this an extremely inefficient method of detection, unless multi-object spectrographs capable of precision radial velocity measurements become available, so that many stars can be monitored simultaneously.

The systems that offer the best chance of observing tidally-captured planets are globular clusters, because they have the highest frequencies of captured planets ($`0.0250.1\%`$), and the largest total number of systems. Furthermore, due to the large velocity dispersions, planets are captured into very close orbits. This increases the a priori probability that the planet will transit its parent stars to nearly unity. In §4, I considered in detail frequency and detectability of tidally-captured planets in 47 Tuc. Although no transiting planets were found in this cluster, it may still be worthwhile to search in other clusters. An observing campaign with the Hubble Space Telescope need not be as ambitious as that toward 47 Tuc.

The expected number of transiting, tidally-captured planets in a rich open cluster is only of order unity. Therefore, the prospects of detecting such planets in transit searches toward open clusters (Street et al., 2003; Burke et al., 2003b) are poor, unless many massive planets are formed per star.

Loose clusters are subject to fairly rapid disruption, after which their stars disperse into the general disk population. It is possible that a large fraction of stars in the disk were initially born in such environments. Assuming a typical lifetime of a loose cluster of $`300\mathrm{Myr}`$, then $`0.005\%`$ of stars in the field should have a tidally-captured planet, if all stars in the field were born in associations, and planet formation were ubiquitous. Accounting for a binary fraction of 50%, and assuming a transit probability of $`40\%`$ for planets captured in loose clusters, I estimate that $`100,000`$ stars need to be monitored to have a chance of detecting even one transiting tidally-captured planet. OGLE monitored 5 million stars toward the Galactic bulge over 45 nights to search for transiting planets to disk stars (Udalski et al., 2002). They only analyze a subset of $`52,000`$ of these stars with photometry better than $`1.5\%`$. Similarly, the EXPLORE project has monitored $`350,000`$ stars in the Galactic plane for 11 nights (Mallén-Ornelas et al., 2003). They achieved better than $`1\%`$ photometry for a subset of 37,000 stars, which they visually inspected for transits. In both cases, an insufficient number of stars were searched to place interesting constraints on the fraction of tidally-captured planets. However, because of the large number of expected transits, it should be possible to detect transiting planets even when the scatter is larger than the expected signal. Therefore, the optimal (but time-consuming) approach would be to search all light curves for transiting planets. For such a photon-limited survey, the detection probability scales as $`P^{5/3}`$ (Pepper, Gould & DePoy, 2003), highly favoring extremely close-in tidally-captured planets.

It is important to reiterate that all of these frequency estimates rely on the assumption that the tidal decay timescale is longer than $`1\mathrm{Gyr}`$, and thus that the quadratic prescription for the suppression of eddy viscosity is approximately correct. The only exception is planets captured in loose associations. Such planets are captured sufficiently far from their parent stars that they are stable over $`\mathrm{Gyr}`$ timescales for both linear and quadratic prescriptions. However, for ages much larger than $`1\mathrm{Gyr}`$, the orbits will rapidly decay under the linear prescription (see Fig. 2). The detection of a planet on an extremely close orbit would therefore have important implications for the tidal dissipation theory.

How can tidally-captured planets be distinguished from planets that migrated to such close orbits? The period distribution of radial-velocity planets shows a pile-up at $`P3\mathrm{d}`$, with a significant lack of planets with smaller periods. Tidally-captured planets would be easily identified in such a distribution. Kuchner & Lecar (2002) suggest that this distribution arises from the fact that the centers of protoplanetary disks are evacuated interior to separations corresponding to periods of $`6\mathrm{d}`$, and that planetary migration thus halts when the planet’s outer Lindlbad resonance reaches the inner disk edge. Although this seems to provide a compelling explanation for the observed period distribution of radial velocity planets, the recent detection of a transiting planet at $`P=1.2\mathrm{days}`$ (Konacki et al., 2003) generally suggests that this picture cannot be universal (Sasselov, 2003). Indeed, Trilling et al. (1998) suggest another mechanism for halting migration: planets migrate until they reach their Roche limit, at which point they lose enough mass to the parent star to balance inward migration. If disk dispersal occurs in a sufficiently short time, the planets will be left intact in short-period orbits. Such planets would generally have similar periods as tidally-captured planets. Therefore, additional diagnostics are needed.

One possible method of distinguishing tidally-captured planets from migrating planets is a radius measurement. Both known transiting planets have radii significantly larger than that of Jupiter. It has been suggested that these large radii are due to the fact that these planets migrated on a sufficiently short timescale that their gravitational contraction was retarded, leaving them in a permanently inflated state (Burrows et al., 2000). If this process is universal, then we can expect migrated planets to have significantly larger radii than tidally-captured planets, which could not have migrated far from their birthplaces. This assumes that there are no processes that can inflate planets once they have been tidally captured. In principle, variable tidal forces from the star on a planet in an eccentric orbit can deposit enough energy to inflate the planet (Gu, Lin, & Bodenheimer, 2003), but to have a significant effect, there must be some external force continuously pumping the eccentricity, since the circularization timescale for extremely close-in planets is very short.

## 6. Summary

Planets which form around stars born in dense stellar environments are subject to dynamical perturbations from other stars in the system. These perturbations will generally serve to strip outer planets from their parent stars, leading to a significant population of freely-floating planets in the system. Further interactions with stars in the system drive this planetary population toward equipartition, thus raising the velocity dispersion of the planets, ultimately resulting in escape from the system. However, some planets will survive long enough to be pass sufficiently close to another star to be tidally captured.

The stability of tidally-captured planets against orbital decay depends critically on the physics of viscous dissipation of the tide in the convective envelope of the star. This processes is extremely uncertain, and different prescriptions lead to decay timescales that differ by four orders of magnitude. For quadratic suppression of eddy viscosity, planets on tidally-captured orbits will generally be stable for $`\mathrm{Gyr}`$, whereas for linear suppression, planets will decay very quickly.

The frequency of tidally-captured planets depends on the competing effects of dynamical stripping, evaporation, capture, and orbital decay. The timescales for these effects depend in turn on the properties of the stellar systems, particularly the velocity dispersion $`\sigma `$ and the number density $`\nu `$. Under the assumption that all stars form four Jupiter-mass planets with a uniform logarithmic distribution in semi-major axis between $`550\mathrm{A}\mathrm{U}`$, and assuming quadratic suppression of eddy viscosity, I have estimated the frequency of tidally-captured planets in a stellar system as a function of the age, velocity dispersion, and number density of the system.

For loose associations with $`\sigma 0.6\mathrm{km}\mathrm{s}^1`$ and $`\nu 10^2\mathrm{pc}^3`$, the frequency of tidally-captured planets after a typical cluster lifetime of $`10^8`$ years is $`0.005\%`$. If most stars formed in such systems, this is roughly the expected frequency of extremely close-in planetary companions to stars in the field. These planets may be found by deep, wide-angle searches for transiting planets around disk stars (Udalski et al., 2002; Mallén-Ornelas et al., 2003). Detection of such a population would imply both that most stars were formed in dense stellar environments, and that the timescale for orbital decay is $`\mathrm{Gyr}`$.

For rich open clusters ($`\sigma 1.5\mathrm{km}\mathrm{s}^1`$ and $`\nu 10^3\mathrm{pc}^3`$), the frequency of tidally-captured planets is $`0.03\%`$. Unfortunately, the prospects for detecting such a population of planets in a given cluster are poor, as $`1`$ transiting planet per cluster is expected.

Due to the long relaxation times and high encounter rates, the frequency of tidally-captured planets in a typical globular cluster ($`\sigma 10\mathrm{km}\mathrm{s}^1`$ and $`\nu 10^4\mathrm{pc}^3`$) is $`0.1\%`$. Therefore, one expects several hundred transiting planets in each globular cluster.

The relatively high frequency of tidal capture in globular clusters makes them excellent targets to search for such systems. The null result of the search for transiting planets in 47 Tuc can be used to place interesting constraints on the initial frequency of planets and brown-dwarfs (either initially free-floating or bound) relative to stars. I find that $`1`$ planet with mass $`M_{\mathrm{Jup}}`$ originally formed per star in 47 Tuc. Less than 50% of stars formed solar-system analogs. The initial frequency of brown dwarfs relative to stars is $`25\%`$. All of these conclusions are under the assumption of quadratic suppression of eddy viscosity. If viscosity suppression is much less efficient, no interesting constraints can be placed.

Tidally-captured planets will have orbital periods of $`0.10.4\mathrm{d}`$, and can be distinguished from migrating planets from their distribution of orbital periods provided the mechanism that is observed to halt planetary migration at $`P3\mathrm{d}`$ is robust. However, if the recent discovery of a planet at $`P1\mathrm{d}`$ is indeed indicating that there exist multiple mechanisms for halting migration, then one cannot rule out the hypothesis that an observed extremely close-in planet was the result of orbital migration. However, the radii of tidally-captured planets are expected to be similar to their radius at their sites of formation, and thus considerably smaller than planets that migrated quickly to close orbits. Thus a radius measurement might provide a diagnostic for tidally-captured planets.

If the majority of stars is formed in dense stellar environments, then most planetary systems will be subject to a variety of dynamical effects that can alter planetary formation, migration, and survival. Therefore, planetary systems cannot be understood as isolated systems, and the interpretation of the observed properties of extrasolar planets will require the consideration of these effects. This study highlights one example of how such considerations can lead to new observable diagnostics of planetary formation.

I would to thank John Bahcall, Ian Bonnell, Neal Dalal, and Roman Rafikov for comments and discussions. This work was supported by NASA through a Hubble Fellowship grant from the Space Telescope Science Institute, which is operated by the Association of Universities for Research in Astronomy, Inc., under NASA contract NAS5-26555.


format: text

text:
# Spectral engineering with multiple quantum well structures

## Introduction

Band-gap engineering is aimed at creating materials with pre-determined electrical properties. From the point of view of optical or optoelectronic applications, it is also important to be able to design materials with pre-determined optical spectra. This requires a development of spectral engineering to control the interaction between light and matter. The realization that such a possibility exists led to the development of the field of photonic crystals.Yablonovitch ; John  Additional opportunities in controlling the light-matter interaction arise in photonic structures made of materials with internal resonances lying in the spectral region of the photonic band structure.DeychPRE1999 ; SivachenkoPRA2001 ; HuangPRL2003  Light propagates through such materials in the form of polaritons, in which electromagnetic waves are coupled with internal excitations of the materials. By changing properties of the material excitations, one can manipulate the properties of light as well. Combining polaritonic effects with photonic crystal effects, one obtains a greater flexibility in designing optical properties and an opportunity to tune them after the growth. An interesting one-dimensional example of such systems is given by Bragg multiple quantum wells (BMQW).IvchenkoMQW  In these structures, the wavelength of the quantum well (QW) exciton radiation, $`\lambda `$, matches the period of the multiple quantum well structure, $`d`$: $`\lambda =2d`$. As a result, the radiative coupling between quantum wells causes a very significant modification of exciton radiative properties, which are effectively controlled by geometrical parameters of the structure. Such structures, therefore, are good candidates for spectral engineering. In Ref. DefectMQW,  it was found that by replacing a single base structural element of BMQW structure with an element with different properties (a defect), one can significantly alter optical spectra of these structures. It was shown that upon introducing different types of the defect, a great variety of spectral types could be created.DefectMQW ; DeychNanotechnology  However, Ref. DefectMQW,  dealt with ideal structures and it was not clear if these effects can be reproduced in realistic structures suffering from homogeneous and inhomogeneous broadenings, and whose lengths are limited by technological capabilities. The goal of this paper is to show that realistic BMQW structures with defects (DBMQW) can be designed to exhibit reflection spectra with sharp features characterized by high contrast even in the presence of relatively large values of the broadenings. We will also show that the spectra of these structures can be tuned after their growth with the help of the quantum confined Stark effect (QCSE).QCSE  This makes DBMQW structures a potential candidate for spectral engineering with applications for tunable switching and modulating devices.

## Reflection spectra of DBMQW structures

We consider a structure consisting of $`N=2m+1`$ QW-barrier layers. The layers are identical except for one in the the middle, where the quantum well has a different exciton frequency. While our calculations are of a rather general nature, we will have in mind a $`GaAs/AlGaAs`$ system as an example. In this case, such a defect can be produced either by changing the concentration of $`Al`$ALdependence\_book ; ALdependence  in the barriers surrounding the central well, or the width of the well itselfwidthdependence  during growth. While both these methods will also affect the optical width of the defect layers, this effect is negligibly small for the systems under consideration, and we will assume that the exciton frequency is the only parameter differentiating the defect well from the others.

The reflection spectra are calculated using the transfer matrix approach. The inhomogeneous broadening of the QW excitons is taken into account within the framework of the effective medium approximation KavokinGeneral  which was shown to describe the main contribution to the reflection coefficient ErementchoukPRB2003  Within this approach the exciton susceptibility, which determines the reflection and transmission coefficients for a single QW, is replaced with its value averaged over the distribution of the exciton frequencies along the plane of a QW:

$$\chi _{h,d}(\omega )=𝑑\omega _0f_{h,d}(\omega _0)\frac{\mathrm{\Gamma }_0}{\omega _0\omega i\gamma }.$$
(1)
Here $`\mathrm{\Gamma }_0`$ is the effective radiative rate of a single QW, characterizing the strength of the coupling between excitons and electromagnetic field, $`\gamma `$ is the parameter of the nonradiative homogeneous broadening, and $`f_{h,d}`$ are distribution functions of the exciton energies of the host and defect QW’s respectively. The variance of this function, $`\mathrm{\Delta }`$, is interpreted as the parameter of the inhomogeneous broadening. The functions $`f_h`$ and $`f_d`$ differ in their mean values, which are $`\omega _h`$ and $`\omega _d`$ for the host and defect QW respectively. The defect-induced effects are most pronounced if $`|\omega _h\omega _d|\mathrm{\Delta }`$. In this case the inhomogeneous broadening of the host wells is negligible in the vicinity of $`\omega _d`$, and defect-induced modifications of the spectra can be studied with only the inhomogeneous broadening of the defect well taken into account.

If the length of the BMQW structure is not very large (for $`GaAs/AlGaAs`$ it should be less than $`500`$ periodsErementchoukPRB2003 ), the reflection coefficient can be presented in the form

$$r=\frac{i\overline{\mathrm{\Gamma }}}{\omega _h\omega +i(\gamma +\overline{\mathrm{\Gamma }})}\frac{\mathrm{\Omega }_s\mathrm{\Gamma }_0D_d}{i\mathrm{\Gamma }_0\mathrm{\Gamma }_0D_d},$$
(2)
where $`D_d=1/\chi _d`$, $`\overline{\mathrm{\Gamma }}`$ is the radiative width of the pure BMQW structure,IvchenkoMQW

$$\overline{\mathrm{\Gamma }}=\mathrm{\Gamma }_0N/(1i\pi qN),$$
(3)
and $`\mathrm{\Omega }_s=(\omega _d\omega _h)/N`$.

The reflection spectrum is characterized by the presence of a minimum and a maximum, Fig. 1, both of which lie in the vicinity of $`\omega _d`$, but are shifted with respect to it. The position of the minimum is determined mostly by parameter $`\mathrm{\Omega }_s`$, $`\omega _{min}=\omega _d\mathrm{\Omega }_s\gamma ^2/\mathrm{\Omega }_s`$. When $`\mathrm{\Omega }_s>\mathrm{\Delta }`$ the inhomogeneous broadening of the defect well does not affect the value of the reflection at $`\omega _{min}`$. For this to happen, the number of wells must satisfy the condition $`N<N_c`$, where $`N_c27`$ for $`GaAs/AlGaAs`$,GaAs  and $`N_c36`$ for $`CdTe/Zn_{0.13}Cd_{0.87}Te`$.CdTe  In this case, the value of the reflection at the minimum $`R_{min}`$ is determined by the rate of the non-radiative relaxation,

$$R_{min}=|\overline{\mathrm{\Gamma }}|^2\gamma ^2N^4/[(\omega _d\omega _h)^4(N1)^2],$$
(4)
and can be very small, when the latter is small. Meanwhile, $`\omega _{max}`$ lies in the spectral region, where the host system is almost transparent, and the value of the reflection at this frequency, $`R_{max}`$ can be estimated as that of a single defect QW

$$R_0=4\mathrm{\Gamma }_0^2/(\pi \stackrel{~}{\gamma }+2\mathrm{\Gamma }_0)^2.$$
(5)
The highest values of the contrast, defined as the ratio of the maximum and minimum reflections $`\eta =R_{max}/R_{min}`$

$$\eta \left[(\omega _d\omega _h)/N\sqrt{\gamma \stackrel{~}{\gamma }}\right]^4$$
(6)
are obtained when the number of periods in the structure is small. For low temperature values of $`\gamma `$, the contrast can be as large as $`10^4`$. However, these large values of the contrast are accompanied by rather small values of $`R_{max}`$. For switching or modulating applications, it would be useful to have large contrast, and a large maximum reflection. The latter can be improved by considering structures with multiple defect wells. This leads, of course, to an increase in the total number of wells, but as we show, one can achieve a significant increase in $`R_{max}`$ for quite reasonable total length of the structure without compromising the contrast too much. Fig. 2 shows the results of numerical computations of the dependence of $`R_{max}`$ and the contrast upon the number of defects. The structures were constructed of several blocks, each of which is a $`9`$-period long BMQW with a single defect well in the middle.

One can see that, indeed, the spectrum of such multi-defect structures exhibits large $`R_{max}`$ (up to $`0.8`$ for structures no longer than $`80`$ periods), while preserving high values of the contrast (of the order of $`10^4`$).

## Tunability

Applications of DBMQW structures for switching or modulating devices is based on the possibility to change the value of the reflection coefficient at a working frequency $`\omega _w`$ by switching between $`\omega _w=\omega _{max}`$ and $`\omega _w=\omega _{min}`$, using for instance the QCSE in order to change the value of $`\omega _d`$. The structures under consideration also allow for tuning of the working frequency of the device by shifting the entire spectrum of the structure using QCSE in host wells. There are several different ways to implement this idea, but here we only want to demonstrate its feasibility. The main difficulty results from the fact that shifting $`\omega _h`$ will detune the whole system from the Bragg resonance and may destroy the desirable spectral features discussed above. In order to see how the detuning affects the spectrum, we assume for simplicity that $`\omega _h`$ and $`\omega _d`$ change uniformly, and study the reflection spectrum of an off-Bragg structure.

It was shown in Ref. DLPRBSpectrum,  that the small detuning from the Bragg resonance results in opening up a propagating band at the center of the forbidden gap significantly complicating the spectrum. It turns out, however, that as long as $`\omega _{min}`$ and $`\omega _{max}`$ are well separated from $`\omega _h`$, the detuning did not affect the part of the spectrum associated with the defect. Indeed, we show that the reflection spectrum of an off-Bragg structure is described by the same Eq. (2) as that of the Bragg structure. The only modification is the change of the definition of $`\overline{\mathrm{\Gamma }}`$, which now becomes

$$\overline{\mathrm{\Gamma }}=\mathrm{\Gamma }_0N/[1iN\mathrm{sin}\pi (\omega \omega _B)/\omega _B].$$
(7)
Thus, for such shifts of the exciton frequencies, $`\omega _s`$, that satisfy the condition

$$N\mathrm{sin}\left(\pi \omega _s/\omega _B\right)1$$
(8)
the destructive effect of the detuning of the structure away from the Bragg resonance is negligible in the vicinity of $`\omega _d`$. It is important to note that the shift should be small in comparison with the relatively big exciton frequency rather than, for example, with the width of the reflection band. Because of this circumstance, our structures can tolerate as large changes of the exciton frequencies as are possible with QCSE. The result of such a change is simply a uniform shift of the part of the spectrum shown in Fig. (1) by $`\omega _s`$.

Additionally, Eq. (7) demonstrates a stability of the considered spectrum with respect to weak perturbations, such as small mismatch of refraction indices of wells and barriers, different optical widths of the host and defect quantum wells, and others.

## Conclusion

In this paper we considered reflection spectra of one special case of DBMQW structures, namely those in which the defect well differs from the host wells in the value of the exciton energy. We showed that if the frequency of the defect lies at the edge of the host reflection band, the spectrum of such a structure becomes significantly modified: in the vicinity of the defect frequency it becomes non monotonic with a well defined minimum and maximum. The value of the reflection at the minimum, $`R_{min}`$ is determined mostly by the rate of the non-radiative relaxation of excitons, and can be very small at low temperatures. The small value of the reflection leads to a giant contrast, defined as $`R_{max}/R_{min}`$, which can be as large as $`10^4`$. The contrast is one of the figures of merit for structures considered for switching or modulating applications, however, the maximum reflection $`R_{max}`$ in such structures is rather low. We showed that $`R_{max}`$ can be significantly increased for structures with several defects without compromising the value of the contrast.

An additional advantage of the proposed structures is their tunability. We demonstrated that shifting the host and defect exciton energies by several widths of the hosts’ reflection band leads to the uniform shift of the entire spectrum without any significant adverse effects on the spectral region in the vicinity of the defect frequency. This shift can be realized using, for instance, QCSE so that the spectra of the considered structures can be electrically tuned.

The work is supported by AFOSR grant F49620-02-1-0305 and PSC-CUNY grants.


text:
# Correlation energy of a two-dimensional electron gas from static and dynamic exchange-correlation kernels

## I Introduction

Fluids of electronic carriers with essentially two-dimensional (2D) dynamics in semiconductor structures present a rich phenomenology, in which dynamical correlations take increasing importance as the electron density is lowered1 . Many of the electron-electron interaction effects can be understood with the help of the homogeneous electron gas (EG) model2 . A central role in the theory of short-range correlations in the EG is played by the wave number and frequency dependent exchange-correlation kernels or equivalently by the local field factors entering the linear response properties of the model3 . These provide key inputs for some applications of density functional theory4  and for studies of quasi-particle properties such as the effective mass and the effective Landè $`g`$-factor5 . These properties are known from experiment for carriers in semiconductor structures over a wide range of carrier density6 . A great deal of accurate information has come from quantum Monte Carlo (QMC) studies7 , but theoretical understanding continues to draw interest.

The exchange-correlation kernel $`f_{\mathrm{xc}}(q,\omega )`$ entering the dielectric response of the EG determines its correlation energy through the fluctuation-dissipation theorem relating the imaginary part of the inverse dielectric function to van Hove’s dynamic structure factor $`S(q,\omega )`$. Starting from $`S(q,\omega )`$ an integration over frequency is needed to obtain the electron-pair structure factor $`S(q)`$, which in Fourier transform describes the shape of the instantaneous exchange-correlation (Pauli-Coulomb) hole surrounding an average electron located at the origin. From $`S(q)`$ the ground-state energy $`E_g`$ can be calculated by means of further integrations over wave number and over coupling strength (see for instance Ref. 3, ) and the correlation energy $`\epsilon _c`$ is defined, as usual, as the difference between $`E_g`$ and its Hartree-Fock value.

For the 3D EG Lein et al.8  have analyzed in detail this procedure for calculating the correlation energy, with particular attention to the roles of the wave number and frequency dependence of $`f_{\mathrm{xc}}(q,\omega )`$. While the wave number dependence of the kernel is physically related to the spatial shape of the exchange-correlation hole, as already remarked, its frequency dependence reflects the inertia of the hole as the electron at the origin moves through the EG. When the frequency dependence of $`f_{\mathrm{xc}}(q,\omega )`$ is omitted and its wave number dependence is approximated by its leading long-wavelength term, one is accounting for short-range exchange and correlation effects in the EG only through the compressibility sum rule (see for instance Ref. 3, ). This approximation provides the basis for the so-called Adiabatic Local Density Approximation (ALDA), which has often been used in dealing with time-dependent phenomena in inhomogeneous electronic systems in the low-frequency regime4 .

In the present work we extend the analysis of Lein et al.8  to the 2D EG by examining how several known forms of $`f_{\mathrm{xc}}(q,\omega )`$ perform in regard to the calculation of its correlation energy over a range of values for the coupling strength. The standard of comparison is the parametrized form reported for the correlation energy by Rapisarda and Senatore9  from their Diffusion QMC results. The paper is organized as follows. In Sec. II we report the detailed expressions relating $`\epsilon _c`$ to $`f_{\mathrm{xc}}(q,\omega )`$ and in Sec. III we list the various forms of $`f_{\mathrm{xc}}(q,\omega )`$ that we have considered. Section IV presents our main numerical results and discusses them with special emphasis on the usefulness of the ALDA in 2D as compared with 3D. Finally, in an Appendix we briefly comment on the scaling properties of the exchange-correlation kernel for inhomogeneous electron fluids of arbitrary dimensionality.

## II Correlation energy from exchange-correlation kernels

We consider a fluid of electrons moving in a plane and interacting by the $`e^2/r`$ law. The correlation energy (per electron) can be written in terms of the difference between the potential energy of the fluid and the exchange energy, integrated over the coupling-strength parameter $`\lambda `$:

$$\epsilon _c=\frac{1}{2}_0^{e^2}\frac{d\lambda }{\lambda }\frac{d^2𝐪}{(2\pi )^2}v_q^\lambda [S_\lambda (q)S_{\mathrm{HF}}(q)].$$
(1)
Here $`v_q^\lambda =2\pi \lambda /q`$ is the 2D Fourier transform of the Coulomb potential with strength $`\lambda `$, $`S_\lambda (q)`$ is the structure factor of an EG with interactions $`v_q^\lambda `$ and $`S_{\mathrm{HF}}(q)`$ is the Hartree-Fock structure factor,

$$S_{\mathrm{HF}}(q)=\{\begin{array}{c}\frac{2}{\pi }\left[\mathrm{arcsin}(\overline{q})+\overline{q}\sqrt{1\overline{q}^2}\right]\text{for }\overline{q}<1\hfill \\ 1\text{for }\overline{q}>1\hfill \end{array}$$
(2)
where $`\overline{q}=q/(2k_F)`$ with $`k_F=(2\pi n)^{1/2}`$ being the Fermi wave number determined by the 2D electron density $`n`$. The fluctuation-dissipation theorem relates the structure factor to the density-density response function $`\chi _{\rho \rho }^\lambda (q,iu)`$ calculated on the imaginary frequency axis,

$$S_\lambda (q)=\frac{\mathrm{}}{\pi n}_0^{\mathrm{}}𝑑u\chi _{\rho \rho }^\lambda (q,iu).$$
(3)
The response function can in turn be expressed through the frequency-dependent exchange-correlation kernel $`f_{\mathrm{xc}}(q,\omega )`$ on the imaginary frequency axis,

$$\chi _{\rho \rho }^\lambda (q,iu)=\frac{\chi _0(q,iu)}{1[v_q^\lambda +f_{\mathrm{xc}}^\lambda (q,iu)]\chi _0(q,iu)}$$
(4)
where

$`\chi _0(q,iu)={\displaystyle \frac{m}{\pi \mathrm{}^2}}\left[1{\displaystyle \frac{1}{\sqrt{2}\overline{q}}}\sqrt{f(\overline{q},\overline{u})+\sqrt{f^2(\overline{q},\overline{u})+4\overline{q}^2\overline{u}^2}}\right]`$ (5)
is the density-density response function of the 2D ideal Fermi gas10 . Here, $`\overline{u}=mu/(\mathrm{}k_Fq)`$ and $`f(\overline{q},\overline{u})=\overline{q}^2\overline{u}^21`$.

From Eqs. (3) and (4) it is possible to rewrite Eq. (1) in the form

$$\epsilon _c=\frac{\mathrm{}}{4\pi ^2n}_0^{\mathrm{}}q𝑑q_0^{\mathrm{}}𝑑u_0^{e^2}\frac{d\lambda }{\lambda }v_q^\lambda \frac{[\chi _0(q,iu)]^2\left[v_q^\lambda +f_{\mathrm{xc}}^\lambda (q,iu)\right]}{1[v_q^\lambda +f_{\mathrm{xc}}^\lambda (q,iu)]\chi _0(q,iu)}.$$
(6)
In the paramagnetic case of present interest $`\epsilon _c`$ depends only on $`n`$, or equivalently on the dimensionless density parameter $`r_s`$ defined by $`r_sa_B=(\pi n)^{1/2}`$, $`a_B`$ being the Bohr radius.

In fact, the full dependence of $`f_{\mathrm{xc}}^\lambda (q,iu)`$ on $`\lambda `$ is not needed for calculating $`\epsilon _c`$ at any given $`r_s`$. It suffices to know the $`r_s`$ dependence of the exchange-correlation kernel at full coupling strength, i.e. $`f_{\mathrm{xc}}(q,iu)f_{\mathrm{xc}}^{\lambda =e^2}(q,iu)`$. That is, the correlation energy can be calculated from

$$\epsilon _c(r_s)=\left[\frac{\sqrt{2}}{r_s^2}_0^{r_s}\gamma _c(r_s^{})𝑑r_s^{}\right]\mathrm{Ryd}$$
(7)
where

$$\gamma _c(r_s)=\frac{1}{k_F}_0^{\mathrm{}}𝑑q[S(q)S_{\mathrm{HF}}(q)]=\frac{\mathrm{}}{\pi nk_F}_0^{\mathrm{}}𝑑q_0^{\mathrm{}}𝑑u\frac{[\chi _0(q,iu)]^2\left[v_q+f_{\mathrm{xc}}(q,iu)\right]}{1[v_q+f_{\mathrm{xc}}(q,iu)]\chi _0(q,iu)}.$$
(8)
Several known approximations to $`f_{\mathrm{xc}}(q,iu)`$ as needed for this calculation will be recalled in the next Section.

## III Approximate kernels for the homogeneous electron gas

We have tested five main forms of the kernel $`f_{\mathrm{xc}}(q,iu)`$. In three of these (designated by the acronyms ALDA, DPGT, and STLS) the frequency dependence of the kernel is omitted. A dynamical kernel has been used in the approximations designated by QSTLS and AKA. Some details are as follows.

(i) ALDA: for the homogeneous electron gas the adiabatic local density approximation omits both the wave number and the frequency dependence of the kernel4 . It amounts to taking into account the compressibility sum rule by setting

$$f_{\mathrm{xc}}^{\mathrm{ALDA}}(q,iu)=\frac{1}{n^2\kappa _0}\left(1\frac{\kappa _0}{\kappa }\right)$$
(9)
where $`\kappa `$ is the compressibility of the EG and $`\kappa _0`$ that of the ideal Fermi gas, given in 2D by $`\kappa _0=\pi r_s^4/2`$ in units of $`a_B^2/\mathrm{Ryd}`$. The ratio $`\kappa _0/\kappa `$ is related to the correlation energy (in Ryd units) by

$$\frac{\kappa _0}{\kappa }=1\frac{\sqrt{2}}{\pi }r_s+\frac{r_s^4}{8}\left[\frac{d^2ϵ_c}{dr_s^2}\frac{1}{r_s}\frac{dϵ_c}{dr_s}\right].$$
(10)
We evaluate this input from the parametrized form of $`\epsilon _c(r_s)`$ of Rapisarda and Senatore9 .

(ii) DPGT: an analytical expression for the static kernel $`f_{\mathrm{xc}}(q,0)`$ as a function of $`r_s`$ in the range $`0r_s10`$ has been obtained by Davoudi et al.11  by fitting the available Diffusion Monte Carlo data for the local field factor in the 2D EG12 . Their expression embodies the compressibility sum rule as well as the asymptotic high-$`q`$ behavior of the static kernel. By setting

$$f_{\mathrm{xc}}^{\mathrm{DPGT}}(q,iu)=f_{\mathrm{xc}}^{\mathrm{QMC}}(q,0)$$
(11)
we are evidently omitting the frequency dependence of the kernel but taking into account its ”exact” dependence on wave number at low frequency.

(iii) STLS: the Singwi-Tosi-Land-Sjölander approximation13  omits the frequency dependence of the kernel by relating it self-consistently to the structure factor $`S(q)`$ through

$$f_{\mathrm{xc}}^{\mathrm{STLS}}(q,iu)=\frac{1}{n}\frac{d^2𝐤}{(2\pi )^2}\frac{𝐪𝐤}{q^2}v_k\left[S(|𝐪𝐤|)1\right].$$
(12)
This expression was justified by an analysis of kinetic equations in the presence of a time- dependent weak external potential and has been used rather widely in the literature.

(iv) QSTLS: a frequency dependence of the local field factor in an STLS-type self-consistent theory has been included by later authors. We adopt in particular the expression14

$$f_{\mathrm{xc}}^{\mathrm{QSTLS}}(q,iu)=\frac{1}{n}\frac{d^2𝐤}{(2\pi )^2}\frac{\chi _0(𝐪,𝐤;iu)}{\chi _0(q,iu)}v_k\left[S(|𝐪𝐤|)1\right]$$
(13)
where the function $`\chi _0(𝐪,𝐤;iu)`$ is given by15

$$\chi _0(𝐪,𝐤;iu)=\frac{m}{\pi \mathrm{}^2}c(1r).$$
(14)
Here $`c=𝐪𝐤/q^2`$ and

$$r^2=\frac{y^2+z^2c^2}{2c^2}+\frac{\left[(y^2+z^2c^2)^2+4c^2z^2\right]^{1/2}}{2c^2}$$
(15)
with $`y=2k_F/q`$ and $`z=2mu/\mathrm{}q^2`$.

(v) AKA: Atwal et al.16  have recently proposed an analytic expression for the dynamical local field factor which is based on many-body perturbation theory and on known sum rules, in the same spirit as in the work of Richardson and Ashcroft17  on the 3D EG. With the notations of Ref. 16, ,

$$f_{\mathrm{xc}}^{\mathrm{AKA}}(q,iu)=v_q\left[G_s(q,iu)+G_n(q,iu)\right]$$
(16)
where

$$G_s(q,iu)=\frac{a_s(u)q+b_s(u)q^7}{1+c_s(u)q+d_s(u)q^7}$$
(17)
and

$$G_n(q,iu)=\frac{a_n(u)q+b_n(u)q^7}{1+c_n(u)q+d_n(u)q^6}.$$
(18)
Here $`a_{s,n}(u)`$, $`b_{s,n}(u)`$, $`c_{s,n}(u)`$ and $`d_{s,n}(u)`$ are functions of frequency and depend on two input parameters, which are the contact value $`g(0)`$ of the pair distribution function and the correlation energy. In our calculations we have replaced Eqs. (41) and (52) of Ref. 16,  for $`d_{s,n}(u)`$, which contain typographical errors, by the following expressions communicated to us by Dr. G. S. Atwal:

$$d_s(\omega )=\frac{\zeta _s\lambda _s^0}{6[\zeta _s\lambda _s^01+g(0)](1+\omega ^4)}+\frac{\lambda _s^{\mathrm{}}\omega }{[1g(0)](1+\omega ^4)},$$
(19)
$$d_n(\omega )=\frac{\zeta _n\lambda _n^0}{(5\zeta _n\lambda _n^0+2\lambda _n^{\mathrm{}})(1+\omega ^4)}+\frac{\gamma _n\omega ^{9/2}}{[5\gamma _n\omega ^{1/2}+2\lambda _n^{\mathrm{}}(1+0.38\omega )^{1/2}](1+0.38\omega )^3(1+\omega ^4)}.$$
(20)
For expressions of the other functions and of the various parameters the reader is referred to the original paper.

In determining the input for this form of the kernel we have again used the form of Rapisarda and Senatore9  for $`\epsilon _c(r_s)`$ and followed Atwal et al.16  in taking the values of $`g(0)`$ from an interpolation formula proposed by Polini et al.18 . We have also checked simple variants of this scheme, taking $`g(0)`$ from calculations in the ladder approximation19  or dropping the frequency dependence of $`f_{\mathrm{xc}}^{\mathrm{AKA}}(q,iu)`$.

## IV Numerical results and discussion

We turn to a presentation of our numerical results, which are collected in Figures 1-4. In Figure 1 (top panel) we show the difference between the correlation energy calculated by means of Eqs. (7) and (8) and the ”exact” correlation energy $`\epsilon _c^{\mathrm{QMC}}`$ from the parametrization of Ref. 9, , for electron densities in the range $`0.5r_s10`$. In 2D this range extends into the intermediate-to-strong coupling regime. In the bottom panel of Figure 1 we report the calculated correlation energies in the weak-coupling regime $`0r_s0.8`$, to show that none of the theories that we have examined satisfies the exactly known limiting expression20  for $`r_s0`$,

$$\epsilon _c(r_s)0.385\frac{2\sqrt{2}}{3\pi }(103\pi )r_s\mathrm{ln}r_s+o(r_s).$$
(21)
This panel also reports the correlation energy given by the Random Phase Approximation (RPA), which corresponds to setting the exchange-correlation kernel to zero and evidently is badly in error even at low $`r_s`$21 .

Turning to the results shown in the top panel of Figure 1, a pleasant surprise in 2D is that ALDA does not overestimate badly the correlation energy as was shown by Lein et al.8  to be the case in 3D. In fact, in the range $`2<r_s<10`$ the ALDA results are seen to be quite close to the DPGT results, in which the full nonlocality of the static kernel in space has been included. The reason for this can be inferred from the work of Davoudi et al.22 , who showed that in 2D satisfying the compressibility sum rule reproduces quite accurately the $`q`$-dependence of the static kernel almost up to $`q2k_F`$. Following Lein et al.8  we have also examined frequency-dependent forms of the long-wavelength kernel and in particular that proposed by Qian and Vignale23 . We have found that this overestimates the correlation energy by as much as $`0.1`$ Ryd at $`r_s=3`$. It appears, therefore, that nonlocality of the kernel in space should be included if its nonlocality in time is accounted for in the calculation of the correlation energy.

This is, of course, what is being done in the AKA calculation. The average deviation of its results in the top panel of Figure 1 is approximately $`0.02`$ Ryd, which is somewhat better than in ALDA or in DPGT. Taking $`g(0)`$ in the AKA kernel from the ladder-diagram calculation19  raises this deviation to about $`0.03`$ Ryd and, most importantly, dropping the frequency dependence of the AKA kernel brings it to about $`0.04`$ Ryd. The summary conclusion of all these calculations is, therefore, that as far as the correlation energy in 2D is concerned one does already fairly well by taking a constant exchange-correlation kernel adjusted to the compressibility sum rule and that nonlocality of the kernel in both space and time needs including in order to obtain some improvement in the results.

Figure 1 in its top panel further shows that STLS gives a very good estimate of the correlation energy, with a deviation of about $`0.004`$ Ryd on average. In this case inclusion of the frequency dependence of the kernel by means of the QSTLS recipe has only a minor effect, bringing the average deviation down to about $`0.0035`$. Consistency of the exchange-correlation kernel with the electron-pair structure thus appears to be important, as far as the calculation of the correlation energy is concerned. This fact may give a useful suggestion for further improvement of the more sophisticated kernels that have been proposed in the more recent work. We may also remark that Dobson et al.24  have recently extended the STLS scheme to inhomogeneous electronic systems and shown that it gives good results for the correlation energy of jellium slabs of finite thickness, lying within $`3\%`$ of the Diffusion Monte Carlo results of Acioli and Ceperley25 .

We turn now to Figures 2-4, in which following again Lein et al.8  we provide an analysis of the correlation energy at $`r_s=1`$ into contributions from correlations between density fluctuations of different wave vectors and different imaginary frequencies. Equations (7) and (8) naturally define a wave-vector analysis $`\epsilon _c(q)`$ if we write

$$\epsilon _c=_0^{\mathrm{}}\epsilon _c(q)d(q/k_F)$$
(22)
and an imaginary-frequency analysis $`\epsilon _c(u)`$ if we write

$$\epsilon _c=_0^{\mathrm{}}\epsilon _c(u)d(2mu/\mathrm{}k_F^2).$$
(23)
Gori-Giorgi et al.26  have obtained an ”exact” wave-vector analysis through the Fourier transform of the coupling-averaged correlation-hole density and built an analytic model for it by combining exactly known limiting behaviors with Diffusion Monte Carlo data, as in the work of Perdew and Wang27  on the 3D EG.

Figure 2 compares the results for $`\epsilon _c(q)`$ given by the theories of present interest with the QMC ones that have kindly been communicated to us by Dr. P. Gori-Giorgi26 . While all theories perform similarly well in this test, it is pleasant to notice the reasonably good behavior shown by the ALDA. Finally, the imaginary-frequency analysis of the various theories is reported in Figures 3 and 4, in the low-$`u`$ and large-$`u`$ regimes respectively. As in the 3D case8 , in all theories $`\epsilon _c(u)`$ starts with a finite negative value at $`u=0`$ and vanishes at large $`u`$. In ALDA $`\epsilon _c(u)`$ becomes positive at $`u7\mathrm{}k_F^2/(2m)`$ and ultimately vanishes from above.

To conclude our work, in Appendix A we briefly comment on how the coupling-constant dependence of the frequency-dependent exchange-correlation kernel in an inhomogeneous electronic system of arbitrary dimensionality may be found from the knowledge of its density dependence at full coupling strength.

###### Acknowledgements.

This work was partially supported by MIUR through the PRIN2001 Initiative. One of us (M. P.) acknowledges partial support from Fondazione Angelo Della Riccia (Firenze) and from the National Science Foundation under Grant DMR 0115974. M. P. T. wishes to thank the Condensed Matter Group of the Abdus Salam International Center for Theoretical Physics for their hospitality during the preparation of the manuscript. We are grateful to Dr. G. S. Atwal for communicating to us Eqs. (19) and (20) and to Dr. P. Gori-Giorgi for sending to us prior to publication the ”exact” results shown in Figure 2.

## APPENDIX A: SCALING OF THE INHOMOGENEOUS EXCHANGE-CORRELATION KERNEL

From time-dependent density functional theory and dynamical scaling Lein et al.8  have derived a coordinate-scaling relation in dimensionality D=3 for the coupling-constant dependence of the exchange-correlation kernel $`f_{\mathrm{xc}}[n](𝐫,𝐫^{};\omega )`$ of an inhomogeneous many-electron system.

It is easy to generalize their proof to an inhomogeneous system in arbitrary space dimensionality D: namely, it is possible to prove that

$$f_{\mathrm{xc}}^\lambda [n](𝐫,𝐫^{};\omega )=\lambda ^2f_{\mathrm{xc}}[n^{}](\lambda 𝐫,\lambda 𝐫^{};\omega /\lambda ^2)$$
(24)
where $`n^{}(𝐫,t)=\lambda ^\mathrm{D}n(𝐫/\lambda ,t/\lambda ^2)`$. In fact, the dimensionality of the system enters the proof only insofar as integrations of the time-dependent density $`n(𝐫,t)`$ over space are involved.


text:
# Ultra-low dissipation Josephson transistor

## Abstract

A superconductor-normal metal-superconductor (SNS) transistor based on superconducting microcoolers is presented. The proposed 4-terminal device consists of a long SNS Josephson junction whose N region is in addition symmetrically connected to superconducting reservoirs through tunnel barriers (I). Biasing the SINIS line allows to modify the quasiparticle temperature in the weak link, thus controlling the Josephson current. We show that, in suitable voltage and temperature regimes, large supercurrent enhancements can be achieved with respect to equilibrium, due to electron “cooling” generated by the control voltage. The extremely low power dissipation intrinsic to the structure makes this device relevant for a number of electronic applications.

In the last few years thermoelectric cooling has gained an increased attention in solid-state physics both from the basic research point of view and as a tool for cryogenic applications, in particular for cooling space-based infrared detectors cryo ; luukanen . The advantage of electronic microrefrigerators stems from their intrinsic compactness, little power dissipation and ease of on-chip integration. To date several refrigeration schemes at cryogenic temperatures have been proposed QD ; nahum ; FS ; paterson ; kap . In particular, cooling schemes based on superconductor-insulator-normal metal (SIN) tunnel junctions proved extremely effective in refrigerating both electrons nahum ; leivo  and phonons manninen . In such systems the cooling principle is based on the superconducting energy gap ($`\mathrm{\Delta }`$) that allows to remove quasiparticles with higher energies from the normal electrode more effectively than those with lower energy averin .

In this letter we propose a device that exploits thermoelectric cooling to obtain transistor effect by means of a symmetric SIN double-microcooler integrated with a SNS weak link (see Fig.1). In existing *controllable* Josephson junctions critical supercurrent modulation is achieved by driving the quasiparticle distribution out of equilibrium vanWees ; wilhelm ; yip  by dissipative current injection into the N region of the junction morpurgo ; baselmans ; shaik ; birge . This leads to controlled supercurrent suppression. Our scheme exploits the exponential temperature dependence of the critical supercurrent in long SNS weak links zaikin ; wilhelm2 . We show that thanks to electric cooling, critical supercurrent values can be strongly increased from their equilibrium magnitude yielding a substantial current gain in the final device. Importantly, thanks to its unique working principle our device leads to control currents that are orders of magnitude lower than in all-normal control lines. This property makes this device of particular interest for extremely-low dissipation cryogenic applications.

The 4-terminal device (sketched in Fig. 1) consists of a N region laterally connected (in the $`y`$-direction) to two superconductors (SJ) through perfectly transmitting interfaces interface , thus realizing a SNS junction of length $`L_J`$ and width $`W_J`$. We consider a mesoscopic diffusive SNS weak link (i.e., $`\mathrm{}_m<L_J<\mathrm{}_\phi `$, where $`\mathrm{}_m`$ and $`\mathrm{}_\phi `$ are the elastic mean free path and the single-particle dephasing length, respectively) and, for reasons that will be clear in the following, we assume that the junction is *long* (i.e., $`\mathrm{\Delta }_J\mathrm{}𝒟/L_J^2=E_{Th}`$, where $`\mathrm{\Delta }_J`$ is the SJ energy gap, $`𝒟`$ is the N region diffusion constant and $`E_{Th}`$ is the Thouless energy). Along the $`x`$-direction the N region is connected to two additional superconducting reservoirs (SR) through identical tunnel barriers (I) of resistance $`_T`$, thus realizing a SINIS symmetric double-refrigerator (cooling control line) leivo ; nahum  of length $`L_R`$ and width $`W_R`$. Resistance $`_T`$ is taken to be much larger than the N region resistance in the $`x`$-direction so that no voltage drop will occur in the N region, thus optimizing the cooling power of the SINIS device averin . We also assume, following Ref. tero , that $`W_RL_J`$ and $`W_JL_R`$. Moreover, we suppose that SJ are kept at the same potential, chosen to be zero, while SR are kept at equal and opposite voltages $`V=\pm V_R/2`$.

We are interested in the regime where the quasiparticle distribution retains a local thermal equilibrium manifesting itself as a Fermi-Dirac function with effective electron temperature $`T_e`$ differing in general from the lattice temperature $`T_l`$ baselmansthesis . This regime can be attained by assuming a large enough $`L_R`$ (i.e., $`\mathrm{}_{ee}<L_R<\mathrm{}_{eph}`$, where $`\mathrm{}_{ee}`$ and $`\mathrm{}_{eph}`$ are the electron-electron and electron-phonon scattering lengths, respectively) so that injected quasiparticles may thermalize through electron-electron interaction. For this to be true we need to focus on lattice temperatures $`T_l<1`$ K, where electron-electron scattering is the main relaxation mechanism. For $`L_R\mathrm{}_\phi ,\mathrm{}_{ee}`$, the cooling control line can be considered as composed by two SIN junctions incoherently added in series. Under bias voltage $`V_R`$ the power $`𝒫_N`$ transferred from the N region to the superconductors SR can be calculated along the lines of Ref. averin . The main results of this analysis are that $`𝒫_N`$ is a symmetric function of $`V_R`$, with a maximum slightly below $`V_R=\pm 2\mathrm{\Delta }_R/e`$. As a function of the lattice temperature, $`𝒫_N`$ shows a maximum of $`0.12\mathrm{\Delta }_R^2/e^2_T`$ at $`T_l0.3\mathrm{\Delta }_R/k_B`$, quickly decreasing at both lower and higher temperatures. As a consequence, both small $`\mathrm{\Delta }_R`$ and low $`_T`$ are desirable in order to achieve high cooling power at sub-Kelvin temperatures. The final $`T_e`$ in the weak link is determined by considering those mechanisms that can transfer energy into the N region. In our case the main contribution is related to electron-phonon scattering given by $`𝒫_{el}=\mathrm{\Sigma }𝒱(T_e^5T_l^5)`$ urbina , where $`𝒱`$ is the N region volume and $`\mathrm{\Sigma }`$ is a material-dependent parameter. The effective temperature $`T_e(V_R)`$ is then determined solving the energy-balance equation $`𝒫_N(V_R,T_e,T_l)+𝒫_{el}(T_e,T_l)=0`$. In writing the above expression we set the quasiparticle temperature in the superconducting electrodes equal to $`T_l`$. Although this is an idealized assumption we expect it, however, to be a reasonable approximation in real devices exploiting quasiparticle traps anghel .

The supercurrent in the diffusive SNS junction ($`I_J`$) can be calculated from the equation yip ; wilhelm ; tero :

$$I_J=\frac{\sigma _N𝒜}{eL_J}_0^{\mathrm{}}f_A(\epsilon ,T_e)j_S(\epsilon ,\varphi ,T_l)𝑑\epsilon ,$$
(1)
where $`\sigma _N`$ is the conductivity of the weak link, $`𝒜`$ its cross section and $`\varphi `$ is the phase difference between the superconductors. According to (1) $`I_J`$ is expressed in terms of the spectral supercurrent $`j_S(\epsilon ,\varphi ,T_l)`$ weighted by the antisymmetric part of the quasiparticle distribution function in the weak link $`f_A(\epsilon ,T_e)`$ (note that the symmetric part does not contribute since $`j_S`$ is an odd function of energy). In the present case of quasi-equilibrium the latter reduces to $`f_A(\epsilon ,T_e)=\mathrm{tanh}[\epsilon /2k_BT_e]`$. It is worth noticing that $`j_S`$ depends on the lattice temperature $`T_l`$ through $`\mathrm{\Delta }_J`$, while $`f_A`$ depends on the electron effective temperature $`T_e`$ in the weak link. The working principle of our controllable junction is described in a simple way in Fig. 2(a, b), where both $`j_S`$ (dash-dotted line) and $`f_A`$ (solid line) are plotted on the same graph as a function of normalized energy. Hatched areas represent the integral of their product, composed of a positive and a negative contribution. In Fig. 2(a) $`T_e`$ is set equal to $`T_l`$, obtained with $`V_R=0`$, while in Fig. 2(b) $`T_e`$ is lowered below $`T_l`$ by biasing the control line, therefore sharpening the quasiparticle distribution. As a result, the supercurrent is enhanced since the positive area is increased expression . This is demonstrated in the following by numerically solving the energy-balance equation for $`T_e(V_R)`$ and inserting it in the expression for $`I_J(V_R)`$ in (1).

In order to simulate a realistic device we choose Al for the SR reservoirs ($`\mathrm{\Delta }_R=180`$ $`\mu `$eV ), Cu for the N region, $`L_R=10`$ $`\mu `$m, $`W_R=0.5`$ $`\mu `$m, $`L_J=1`$ $`\mu `$m and $`W_J=0.4`$ $`\mu `$m. Supposing a 100-nm-thick film, the total volume of N region amounts to $`𝒱=0.52`$ $`\mu `$m3. For metals like Ag, Cu, and Au the chosen value for $`L_R`$ corresponds to the interacting hot-electron regime morpurgo . Furthermore we take $`\mathrm{\Sigma }=2\times 10^9`$ W K-5 $`\mu `$m-3 nahum  and $`_T=100`$ $`\mathrm{\Omega }`$ blamire . For SJ we choose Nb ($`\mathrm{\Delta }_J=1.52`$ meV) and assuming a diffusion constant $`𝒟=0.007`$ m2/s pierre , the Thouless energy is $`E_{Th}4.6`$ $`\mu `$eV, so that the considered junction is in the long limit ($`\mathrm{\Delta }_J/E_{Th}330`$) nbnote .

In Fig. 3 we plot the calculated normalized $`I_J(V_R)`$ characteristic for several starting lattice temperatures ($`_N=L_J/\sigma _N𝒜`$ is the SNS junction normal state resistance). For all chosen $`T_l`$, $`I_J(V_R)`$ increases monotonically for bias up to $`eV_R2\mathrm{\Delta }_{Al}`$, for which the effect of the control line is to decrease the effective electron temperature. In the case of an *all-normal* control line (i.e., without SR) the effect would have been only to enhance the electron temperature ($`T_e>T_l`$) morpurgo  with consequent suppression of the supercurrent. For larger bias, in fact, an efficient suppression of the supercurrent occurs due to electron heating. As a result the device behaves as a fully-tunable superconducting transistor. At zero temperature a critical current of about $`54\mu `$A can be achieved (it can be further enhanced in a slightly shorter weak link, resulting into a larger $`E_T`$ and a smaller $`_N`$ tero ). Remarkably, the current gain $`𝒢=\delta I_J/I_{Rmax}`$, i.e., the ratio between the maximum supercurrent enhancement at a given $`T_l`$ and the injected quasiparticle current needed to produce such enhancement, obtains values as large as $`𝒢90`$ in the examined $`T_l`$ range. Note that $`𝒢`$ is dependent on the choice of the structure parameters.

We can now clarify the reason why *long* junctions are more appropriate for the device to work. Since in long junctions tero  $`E_{Th}`$ sets the temperature scale for the supercurrent, for sufficiently low lattice temperatures $`I_J`$ is independent of $`T_l`$ and only depends on $`T_e`$ as displayed in the inset of Fig. 3: by starting from a given temperature $`T_e`$ at $`V_R=0`$, an increase of $`V_R`$ lowers $`T_e`$ thus enhancing $`I_J`$. This shows that it is advisable to choose a temperature window where the slope of the $`I_J(T_e)`$ characteristic is large in the range of a high cooling power from the control line, in order to achieve large supercurrent enhancements (of the order of several hundreds percent) with small temperature variations. This optimized temperature range is marked by the hatched area and is well within the cooling performance of the chosen refrigerators. Note that in the *short* junction limit ($`\mathrm{\Delta }_JE_{Th}`$), it is the energy gap $`\mathrm{\Delta }_J(T)`$ which sets the temperature dependence of the supercurrent, thus implying a much reduced effect from the cooling line.

We would like, furthermore, to comment on the impact of power dissipation in this device for practical applications. For fixed geometry a straightforward comparison shows that, in the present case, control currents $`I_R`$ corresponding to bias larger than $`2\mathrm{\Delta }_R/e`$ are at least two orders of magnitude lower ($`I_R10^6`$ A) corrente  than in the existing *all-normal* controllable Josephson junctions morpurgo ; baselmans ; shaik ; birge . For bias corresponding to supercurrent enhancement, on the other side, control currents are at least four orders of magnitude lower, due to the dependence of tunneling subgap current from barrier resistance. These considerations suggest that such a device can be well suited for high-density and ultra-low dissipation cryogenic applications.

In conclusion, a promising scheme to obtain transistor effect in a Josephson weak link was presented. A realistic Nb/Cu/Nb weak link in combination with Al/Al2O3/Cu refrigerators was analyzed and we demonstrated that a supercurrent enhancement of several hundreds percent can be achieved. Moreover, the drastic reduction of control currents intrinsic to the proposed design with respect to existing schemes makes such a structure of interest for applications where extremely low power dissipation is required at cryogenic temperatures.

The authors gratefully acknowledge J. P. Pekola for fruitful discussions. The present work was supported in part by INFM under the PAIS project TIN.