Eccentric Stellar Discs with Strong Density Cusps and Separable Potentials

CYCLE SLIPS IN THE PLL SYNCHRONIZER WITH LIMITER PHASE DETECTOR IN THE PRESENCE OF NOISE ANDCOCHANNEL INTERFERENCEPavel HasanTelecommunications Institute University of Erlangen-Nuremberg Cauerstr. 7, D-91058 Erlangen, GermanyandInstitute of Radio Engineering and Electronics Academy of Sciences of the Czech Republic Chaberska 57, CZ-18251 Prague 8, Czech RepublicAbstract: The paper addresses the cycle slip performance analysis of the phase-locked loop (PLL) synchronizer with limiter phase detector in the presence of noise and cochannel interference. Conditioned on the interference phase, an equivalent loop model of the limiter phase detector is developed under this situation. The conditional cycle slip rate is obtained by adapting a recent formula, and the mean cycle slip rate is reported for the uniformly distributed interference phase. It is shown that the cochannel interference enhances the mean cycle slip rate by reducing the relevant effective loop signal-to-noise ratio.slip statistics of the considered synchronizer is compared to that of a first-order loop with the multiplier phase detector.II. MODELLINGConsider the sinusoidal type limiter phase detector implemented by a multiplier and a bandpass limiter in the signal path [4] as a part of a first-order phase-locked loop.According to the channel model under consideration, the desired signals t ()=2S sin 0t +Ω0t ()(1)Keywordsto be phase tracked by the PLL synchronizer is additively disturbed by AWGN n t () and cochannel interferenceSynchronization, phase-locked loops, cochannel interferenceI. INTRODUCTIONPhase-locked loop (PLL) synchronizers with hard limiters preceding the phase detectors are of great practical interest.As manifestation of their inherently nonlinear nature, cycle slips appear at low loop signal-to-noise ratios. Performance analysis of phase-locked loops in the presence of additive white Gausssian noise (AWGN) [1], [2] shows the cycle slip statistics of the limiter phase detectors superior to that of the conventional multiplier phase detector. Although the presence of a cochannel interferer has become a very probable event, little seems to be known about the PLL synchronizer performance under this situation. Kostic [3]investigated the effects of AWGN and cochannel interference on the characteristics of limiter phase detectors in the open-loop condition. Due to the state-dependent loop noise spectral density of the limiter phase detectors, it is,however, not possible from their characteristics only to deduce the closed-loop PLL performance.I t ()=2J sin 0t +Ω0t +().(2)Here, Ω0 is the open-loop detuning relative to the loop quiescent frequency 0, is the time-invariant interference phase relative to the desired signal with probability density function (pdf) p (), and n t () is a zero-mean stationary narrowband process with two-sided spectral density N 0W/Hz and (one-sided) bandwidth B C >>W L with W L as the interference-free loop noise bandwidth at the operating signal-to-noise ratio (SNR) defined below.Conditioned on the interference phase , the phase detector is modelled [4] by the effective phase detector characteristicg =E sin+([(3)The present paper focuses on the closed-loop performance of limiter phase detectors in the presence of AWGN and cochannel interference. An equivalent loop model of the limiter phase detectors conditioned on the interference phase is developed by adapting the approach originated by Rosenkranz [4]. The model is applied to the sinusoidal type limiter phase detector implemented by a multiplier and a bandpass limiter in the signal path. The conditional cycle slip rate of a first-order loop is evaluated using a recently derived formula [5], and the mean cycle slip rate is given for the uniformly distributed interference phase. The cyclea n d t h e s t a t e -d e p e n d e n t e f f e c t i v e p h a s e n o i s e ′n t ()=sin +(−g )with variance′n 2)=E sin 2+([−g 2).(4)Here, is the loop phase error, is the random phase ofthe sum of the signal, noise and interference, and the conditional expectations E ⋅[] are evaluated using the conditional pdf [6]p =12+a n cos n −()[]n =1∞∑(5)with [8]B N =BC 1+0.098exp −h 121−()[]{}(14)where [7]as the equivalent bandwidth of N t().a n =h exp −h 22()I n −1()2h 22()+I n +1()2h 22()[](6)III. ANALYTICAL RESULTSandFor a given interference phase , the conditional mean cycle slip time T S () is the mean time for the loop to first slip a cycle starting from an arbitrary initial phase error 0,or, equivalently, the mean time between cycle slips. The conditional mean cycle slip time of the considered synchronizer is given by [5]=atanh 2sin1+h 2cos.(7)In (6),h 2=h 121+h 22+2h 2cos()(8)T S ()L (D 0()d0+2∫exp Ux ,()−U ,()−2⌠ dx(15)where h 12=S N 0B C () i s t h e S N R , h 22=J S i s t he interference-to-signal ratio (ISR), and I ⋅() is the modified Bessel function of the first kind and order .Using the above modelling, the effective phase detctor characteristic and the variance of the effective phase noise conditioned on are given bywhere, conditioned on ,W L()=K h 12,h 22,()2(16)g =a 1sin +()(9)andis the effective two-sided loop noise bandwidth,′n 2)=121−2a 12+2a 12−a 2()cos 2+()[]{},(10)()=1W L ()S N 0,0,()(17)respectively. Cycle slips in the first-order phase-locked loop are governed by the stochastic differential equationis the effective loop SNR,D 0()=11+exp 2()()1−a 2a a +b ()(18)=Ω0−K h 12,h 22,()sin+N t ([](11)where t ()=t ()+ is a phase error process with thec y c l e s l i p s t a t i s t i c s e q u a l t o t h a t o f t (), a nd K h 12,h 22,()=a 1K D K 0is the effective loop gain of the linearized loop at the operating SNR and ISR with K D and K 0 as the phase detector gain in the absence of both noise and interference and the gain of the voltage-controlled oscillator, respectively. The conditional loop noiseand U eq x ,()=()1−a 2()a +b a cot x−1ab atan b a cos x(19)N t ()=1(12)withis a within the loop bandwidth white process with two-sided spectral density S N 0,,()=S N f ,,()f =0 given by()=Ω0K h 12,h 22,()(20)S N 0,,()=1N(13)as the conditional normalized effective open-loop detuning,a =1−22a 12+a 2>0 and b =22a 12−a 2()>0.IV. NUMERICAL RESULTSFig.1 shows the conditional cycle slip rate normalized by the loop noise bandwidth W L versus the interference phase for the loop SNR =10dB with the ISR as parameter.Figs. 2 and 3 compare the normalized mean cycle slip rate f o r =0 and 0.5, respectively, of the considered synchronizer to that of a first-order loop with the multiplier phase detector [9], [10] at the same operating loop SNR.The cycle slip statistics of the PLL synchronizer is in accordance with [9] given in terms of the cycle slip rate which appears to be a more appropriate performance measure than the mean cycle slip time under this situation.The conditional cycle slip rateS ()=1T S ()(21)Fig. 2. Normalized mean cycle slip rate versus loop SNR.Limiter phase detector (continuous lines) and multiplier phase detector (dashed lines), =0.obtained by evaluating the double integral in (15)numerically using the Gauss-Legendre quadrature is averaged over p () to yield the mean cycle slip rate E S [].For p ()=12() considered here,E S []=12S ()d −∫.(22)Numerical results reported are given in terms of the loop noise bandwidthW L =K h 12,0,⋅()2,(23)the normalized open-loop detuningFig. 3. Normalized mean cycle slip rate versus loop SNR.Limiter phase detector (continuous lines) and multiplier phase detector (dashed lines), =0.5.=Ω0K h 12,0,⋅()(24)and the loop SNR=1W L S N 0,0,⋅()(25)at the operating SNR in the absence of interference.Fig. 1. Normalized conditional cycle slip rate versus interference phase. =0, =10dB.V. CONCLUSIONSTo author’s knowledge for the first time the closed-loop performance analysis of the limiter phase detector in the presence of AWGN and cochannel interference has been reported. The mean cycle slip rate of the PLL synchronizer with the sinusoidal type limiter phase detector has been evaluated and compared to that of the PLL synchronizer with the multiplier phase detector at the same operatingloop SNR. Assuming the same operating loop noise bandwidth W L and the uniformly distributed interference phase, the results obtained show performance degradation of the former to that of the latter due to the cochannel interference at ISR>-10 dB. A t l o w I S R , t h e P L L synchronizer with the limiter phase detector outperforms that with the multiplier phase detector by several orders of the mean cycle slip rate. The analysis reported provides one with insight into the closed-loop performance of the limiter phase detector under realistic carrier tracking conditions.VII. REFERENCES[1]P. Hasan, "Mean Cycle Slip Time of First-Order Phase-Locked Loop with Limiter Phase Detector" Electron.Lett., vol. 30, no. 5, March 1994, pp. 387-389.[2]P. Hasan, "Approximate Mean Cycle Slip Time of Synchroniser with Triangular Type Phase Detector"Electron. Lett., vol. 30, no. 25, Dec. 1994, pp. 2095-2096.[3]I.M. Kostic, "Effects of Noise and Cochannel Interference on Phase Detector Transfer Characteristic"AEÜ, vol. 43, no. 1, Jan./Feb. 1989, pp. 53-54.The mean cycle slip rate of the PLL synchronizer with the limiter phase detector increases due to the cochannel interference by several orders of magnitude compared to the interference-free value. This is explained by reducing the effective loop SNR () conditioned on as follows. The conditional cycle slip rate S () is a steeply decreasing function of () which is maximum and minimum at =0 a n d =±, respectively. Considering the interference phase distributed uniformly over −,(), the mean cycle slip rate E S []is given largely by the contributions from S () a t close to ±. In fact, the mean cycle slip rate is enhanced by the reduced relevant effective loop SNR due to the cochannel interference at near ±. Morover, at low loop SNR, the loop noise spectral density becomes independent of the phase error, and the limiter phase detector performs like the multiplier phase detector irrespective of the ISR.[4]W. Rosenkranz, "Phase-Locked Loops with Limiter Phase Detectors in the Presence of Noise" IEEE mun., vol. COM-30, no. 10, Oct. 1982, pp. 2297-2304.[5]P. Hasan, "Limiter Effect on Cycle Slip Performance of Synchronization Loops" IEEE Trans. Commun., vol.44, no. 5, May 1996, pp. 562-565.[6] D. Middleton, An Introduction to Statistical Communication Theory, McGraw-Hill, New York: 1960,p. 417.[7]P. Hasan, "Low SNR Approximation to Phase PDF for PSK Signals" Electron. Lett., vol. 30, no. 24, Nov.1994, p. 2008.[8]P. Hasan, "Modelling of Limiter Phase Detector in the Presence of Noise and Interference" Frequenz, vol. 47,no. 5-6, May-June 1993, pp. 129-133.[9]P. Hasan, "Cycle Slip Statistics of a First-Order Phase-Locked Loop in the Presence of Additive Noise and Cochannel Interference" European Trans. Telecommun.,vol. 5, no. 5, Sept.-Oct. 1994, pp. 619-623.VI. ACKNOWLEDGMENT[10]C.Y. Yoon and W.C. Lindsey, "Phase-Locked Loop Performance in the Presence of CW Interference and Additive Noise" IEEE Trans. Commun., vol. COM-30, no.10, Oct. 1982, pp. 2305-2311.T h i s w o r k h a s b e e n s u p p o r t e d b y t h e D e u t s c h e Forschungsgemeinschaft.。

IBM's centenary百年IBMThe test of time时间的考验Which of today’s technology giants might still be standing tall a century after their founding?今天的科技巨头，哪一个也能在它成立后仍续领风骚一百年？010 Leaders - IBM_s centenary.mp3(2.11 MB, 下载次数: 327)Jun 9th 2011 | from the print editionIT IS not, by any means, the world’s oldest company. There are Japanese hotels dating back to the 8th century, German breweries that hail from the 11th and an Italian bank with roots in the 15th. What is unusual about IBM, which celebrates its 100th birthday next week, is that it has been so successful for so long in the fast-moving field of technology. How has it done it?它绝对不是世界上历史最悠久的公司。

日本酒店可追溯到8世纪，德国啤酒厂的诞生于11世纪，而意大利银行起源于15世纪。

下周庆祝它的100岁生日的IBM公司与众不同的是，它长久以来在日新月异的技术领域的已经获得如此成功。

它是怎样做到的呢？IBM’s secret is that it is built around an idea that transcends any particular product or technology. Its strategy is to package technology for use by businesses. At first this meant making punch-card tabulators, but IBM moved on to magnetic-tape systems, mainframes, PCs, and most recently services and consulting. Building a company around an idea, rather than a specific technology, makes it easier to adapt when industry “platform shifts” occur (see article).IBM公司的秘诀在于它秉承着“超越任何具体产品和技术”的理念。

小学下册英语第6单元测验卷[有答案]英语试题一、综合题(本题有50小题，每小题1分，共100分.每小题不选、错误，均不给分)1 I enjoy ______ (与同龄人互动).2 Fish breathe underwater using their __________.3 What is 8 divided by 2?A. 2B. 3C. 4D. 54 What do we breathe?A. WaterB. AirC. FoodD. Light答案: B5 We have ______ (math) class at o'clock.6 A _______ reaction is one that absorbs heat. (吸热)7 The kitten likes to chase a _________. (线)8 It is ________ (sunny) today.9 My dad is very ________.10 My neighbor has a big _______ (我邻居有一只大_______).11 What is the name of the famous wizarding school in Harry Potter?A. BeauxbatonsB. DurmstrangC. HogwartsD. Ilvermorny答案:C12 The ice is very ___. (slippery)13 My uncle is a ______. He enjoys woodworking.14 The whale sings deep in the _________. (海洋)15 The __________ (历史的叙述) tells our stories.16 I have a _______ (surprise) for you.17 My sister is my best ______ (朋友). We play games together and share our ______ (秘密).18 I like to ______ (参加) community projects.19 The capital city of Cuba is __________.20 A saturated solution contains the maximum amount of _____ that can dissolve.21 What is the main ingredient in pizza sauce?A. CheeseB. TomatoC. MeatD. Vegetables答案:B22 Certain plants have developed unique ______ to thrive in their environments. (某些植物发展出独特的特征以适应其环境。

a rXiv:as tr o-ph/2514v22Nov22Eccentricity Evolution of Extrasolar Planets Eccentricity Evolution of Extrasolar Multiple Planetary Systems due to the Depletion of Nascent Protostellar Disks M.Nagasawa 1,2,D.N.C.Lin 1,and S.Ida 3mnagasawa@,lin@,ida@geo.titech.ac.jp ABSTRACT Most extrasolar planets are observed to have eccentricities much larger than those in the solar system.Some of these planets have sibling planets,with com-parable masses,orbiting around the same host stars.In these multiple planetary systems,eccentricity is modulated by the planets’mutual secular interaction as a consequence of angular momentum exchange between them.For mature plan-ets,the eigenfrequencies of this modulation are determined by their mass and semi-major axis ratios.But,prior to the disk depletion,self gravity of the plan-ets’nascent disks dominates the precession eigenfrequencies.We examine here the initial evolution of young planets’eccentricity due to the apsidal libration or circulation induced by both the secular interaction between them and the self gravity of their nascent disks.We show that as the latter effect declines adiabati-cally with disk depletion,the modulation amplitude of the planets’relative phase of periapse is approximately invariant despite the time-asymmetrical exchange of angular momentum between planets.However,as the young planets’orbits pass through a state of secular resonance,their mean eccentricities undergo systematicquantitative changes.For applications,we analyze the eccentricity evolution of planets around Upsilon Andromedae and HD168443during the epoch of proto-stellar disk depletion.We find that the disk depletion can change the planets’eccentricity ratio.However,the relatively large amplitude of the planets’eccen-tricity cannot be excited if all the planets had small initial eccentricities.Subject headings:planetary systems:formation—celestial mechanics—stars:individual(Upsilon Andromeda,HD168443)—planetary systems:protoplane-tary disks—extra solar planets1.INTRODUCTIONMany extrasolar planetary systems have been discovered recently using the radial ve-locity technique(Marcy&Butler2000).The basic assumption is that the spectroscopic variations observed in some targeted stars are due to the Doppler shift associated with the reflex motion of stars with unseen companions.The mass of these companions depends on the poorly known inclination of these systems.But,unless the orbits of these systems are highly inclined(Stepinski&Black2001),the inferred masses of the companions are compa-rable to that of Jupiter.The observed spectra of some targeted stars indicate the presence of multiple companions around them.For example,three companions are found to orbit around Upsilon Andromedae.The long term stability of this system requires their inclination to be sufficiently small such that the masses of these companions are no more than a few Jupiter mass(M J)and much smaller than that of their host star(Laughlin&Adams1999;Rivera &Lissauer2000;Stepinski,Malhotra,&Black2000;Ito&Miyama2001;Lissauer&Rivera 2001).Two companions are also found around HD168443with minimum masses which are an order of magnitude larger than M J.Although the masses of these companions would be a modest fraction of a solar mass(M⊙)if this system is viewed nearly face on,such com-pact hierarchical stellar systems have not be seen before.Thus,we follow the conventional practice to refer these multiple companions as planets.In the limit that the planets’masses are substantially smaller than that of their host star,their mutual secular perturbation induces them to exchange angular momentum while preserving their energy.(Planets in mean motion resonances also exchange energy on com-parable time scales).Around Ups And and HD168443,the planets’orbits are not in mean motion resonances so that their semi-major axes are conserved while their eccentricity and longitude of periapse modulate over some characteristic secular time scale(Murray&Der-mott1999).But,the secular perturbation between the planets has not always been sustained at the present level.During the epoch of their formation,protoplanets are embedded in protostellar disks which have been found around most young stellar objects(cf.Haisch,Lada,&Lada 2001).The mass and temperature distribution of these disks are very similar to those inferred from the minimum mass nebula model for the solar system(Beckwith1999).The self gravity of these disks can induce the orbits of planets formed within them to precess at a rate fasterthan that due their mutual perturbation.Consequently,the rate of angular momentum transfer between the interacting planets is suppressed.In the solar system,the initial contribution of the solar nebula to the total gravita-tional potential dominates the precession frequency of asteroids and comets over that due to the secular perturbation induced on them by the giant planets.The planets also undergo precession induced by the disk gravity and other planets’secular perturbation.In general the precession frequencies of these celestial bodies do not equal each other.But as gas is depleted in the solar nebula,its self gravity weakens.The total precession frequencies of both the asteroids and the planets declines,though at a different rate.When the precession frequency of asteroids with some semi major axes coincides with that of a major planet,they enter into a state of secular resonance.In this resonance,the eccentricity of the asteroids is either excited or damped monotonically,depending on their relative longitude of perias-tron passage with respect to that of the planet.As this secular resonance sweeps across the solar system,large eccentricity may be excited among some small celestial bodies(Ward, Colombo,&Franklin1976;Heppenheimer1980;Nagasawa,Tanaka,&Ida2000;Nagasawa &Ida2000).In this paper,we examine the effects of disk depletion on the eccentricity evolution of the planetary systems around Ups And and HD168443.Through such an investigation,we hope to infer the kinematic properties which these planets are born with and thereby cast constraints on their formation process.In§2,we briefly discuss the precession due to the secular interaction between planets and that due to the self gravity of the protostellar disk. The conditions for secular resonance are discussed.We introduce,in§3,a working model and describe the method we used to analyze the eccentricity evolution.In§4,we present the results of some calculations.We use these numerical results to demonstrate the evolution of these systems with Hamiltonian contour maps.Based on these results,we infer,in§5,some implications on the planets’orbital eccentricity shortly after their formation and while they are still embedded in their nascent disks.Finally,we summarize our results in§6.2.SECULAR INTERACTION IN MULTIPLE EXTRASOLARPLANETARY SYSTEMS2.1.Current Orbital PropertiesIn this paper,we focus our discussion on the planets around Upsilon Andromeda(c and d)and those around HD168443.The decomposition of Upsilon Andromeda’s spectra indicate that there are three planets orbiting around it(Butler et al.1999).Inferred orbitalelements4of the Upsilon Andromeda planets are shown in Table1(a).The semi-major axis,the eccentricity,the longitude of periastron,mass of the planet,and periastron passage time(JD)are denoted by a,e,̟,M,and T peri,respectively.The subscripts b,c,and d representthe values for the individual planets.Throughout this analysis,we assume that these systemsare viewed edge-on and their masses correspond to their minimum values.Provided theirorbits are coplanar,our analysis is independent of the planets’inclination.Our approachbecomes inadequate for the limiting cases of nearly face-on orbits where the masses of thecompanions are comparable to their host stars.Around Ups And,the eccentricity of the innermost planet b is essentially undetectable.The orbit of planet b is likely to be circularized during the main sequence life span of the hoststar Ups And by the tidal dissipation within the planet’s interior as expected for Jupiter-likeextrasolar planets with semi-major axis less than0.05AU(Rasio et al.1996).With its lowmass and small semi-major axis,planet b does not contribute significantly to the dynamicalevolution of the system(Mardling&Lin2003).Although the outer planets(c and d)stillexert secular perturbation on planet b,the cumulative effect on its eccentricity modulationis limited.Under the present configuration,the secular interaction of planet b with planets cand d is weakened by the rapid precession due to the post-Newtonian relativistic correctionin the gravitational potential of the host star(Mardling&Lin2003).Neglecting the secularperturbation due to planet b,the orbital evolution of planets c and d obtained from the directorbital integration of the full equation of motion is shown in Figure1a.In the numericalintegration,the calculation is started with the eccentricity ratio x≡e c/e d and the relativelongitude of periastrons between planets c and d,η≡̟c−̟d,based on their observed values0.66and∼−0.06radian,respectively.We compute the evolution of the system forover2×104years.The equi-Hamiltonian contours(see Appendix D)are also shown.Forplanets c and d around Ups And,the massive planet d has larger eccentricity than smallerplanet c.The longitudes of periastron of planet c and d are always close to each other.Whenthe planets are in librating(closed)track in(x-η)diagram,the planetary system tends to bestable for a long time.The stabilities of this system around Ups And are well investigated(Laughlin&Adams1999;Rivera&Lissauer2000;Stepinski et al.2000;Barnes&Quinn2001;Ito&Miyama2001;Lissauer&Rivera2001;Chiang,Tabachnik,&Tremaine2002;Mardling&Lin2003).Around HD168443,two massive planets are inferred from the radial velocity curves(Marcy et al.2001).The best-fit orbital elements are shown in Table1(b).We also numeri-cally integrate the present-day orbital evolution of the planets b and c around HD168443.Incontrast to the planetary system around Ups And,the lower mass planet b around HD168443 system,has an eccentricity which is more than twice that of the massive planet c.The orbits are integrated over2.5×104years,starting with x=2.65andη∼1.92radian.These planets evolve in circulating(open)track in the(x-η)diagram(Fig.1b)andηevolves from−πtoπover a period of about1.6×104years.This system too is stable despite the large magnitude of eccentricity of planet b.The origin of such large mass and eccentricity of planets in this system has not been addressed previously.2.2.Secular Perturbation between Two PlanetsThe secular interaction induces eccentricity modulation between planets c and d around Ups And(Laughlin&Adams1999;Rivera&Lissauer2000;Ito&Miyama2001;Mardling &Lin2003).For presentation purpose,it is useful to briefly recapitulate the analysis of secular interaction between two planets.It is customary to consider the secular(long-term) evolution of the interacting planets’orbits using a disturbing function(Murray&Dermott 1999).To the lowest order,the modulation of the eccentricity(e c,d)and longitude of periapse (̟c,d)of some planets c and d can be approximated byde c,d=γc,d 1+C e d,cdτhighlight the dominant physical effects.For the numerical calculation of the planets’orbits (see§4),we use the exact equation of planets’motion.Planet-disk interaction also leads to both axisymmetric modification of the gravitational)(Goldreich&Tremaine1980,1982;Lin potential(N c,d)and non axisymmetric torque(Λe,̟c,d&Papaloizou1986a,b,1993).For the latter effect,planets excite waves in the disk which carry angular momentum.The dissipation of these waves,anywhere in the disk,would result in afinite torque(Papaloizou&Lin1984).In principle,the net torque vanishes in the invisid limit.However,their amplitude grows and steepens into nonlinear shocks as the waves propagate away from the location where they are launched,leading to an effective torque(Savonije,Papaloizou,&Lin1994).The evolution of the axisymmetric potential modifies the precession of̟but do not directly influence the eccentricity of the orbits(see eq.[2]).While the non axisymmetric torque may induce monotonic changes in e(Chiang &Murray2002)and̟(see eq.[1]),the magnitude and sign ofΛe c,d depends sensitively on the disk structure.Interaction between the embedded planets with disk gas through corotation resonances damps the eccentricity while that through Lindblad resonances excites the eccentricity(Goldreich&Tremaine1980).For protoplanets with mass less than a few times that of Jupiter and modest eccentricity,gas mayflow in the vicinity of their orbits such that the eccentricity damping effect of the corotation resonances is stronger than the excitation effect of the Lindblad resonances(Goldreich&Tremaine1980).But protoplanets with masses an order of magnitude larger than that of Jupiter may open relatively wide gaps in protostellar disks.In this limit,the protoplanets’corotation resonances may be cleared of disk gas such that their eccentricity may be excited(Artymowicz1993;Papaloizou,Nelson, &Masset2001;Goldreich&Sari2002).For precession,the contribution from the non axisymmetric torque may be weaker than that due to the disk’s axisymmetric contribution to the total gravitational potential in the limit that the torque resulting from the planets’interaction with the interior and exterior regions of the disk are balanced or in low-viscosity and thin disks where planets with relatively low masses can open wide gaps.However,the non axisymmetric torque may lead to significant contributions over time in relatively massive disks.In order to analyze the contribution of each effect,we adopt a piece meal approach.In this paper,we focus our attention on the evolution of planetary orbits due to the changes in the axisymmetric disk potential.For the departure from a point-mass potential,the apsidal motion of planets c and d are included in N c,d(see Appendix A).In thisfirst step,we neglect the effects due to the non axisymmetric torque by settingΛe,̟=0in equations(1)c,dand(2).The extension of our discussion to the limit offinite non axisymmetric torque will be presented in a future contribution.2.3.Precession FrequenciesTo the lowest order,the solution of equations(1)and(2)can be expressed ase c,d exp(i̟c,d)=A c,d exp{i(g1t+β1)}+B c,d exp{i(g2t+β2)},(3) where A c,d and B c,d are the oscillation amplitudes,β1andβ2are the phase angles(e.g., Brouwer&Clemence1961).The individual planet’s longitudes of periastrons precess with two independent eigenfrequencies,g1,2=1t c,g cd≡Cγ1/2dτc 2+4g cd 1/2,(6)which vanishes on the exact center of secular resonance(where g1=g2)only if g c=g d and g cd=0.The former requirement corresponds to a necessary resonance condition A1=0whereA1≡1−γd+∆N(7) is the precession rate induced by axisymmetric component of the perturbed potential(see eq.[9]below).The latter requirement is satisfied if either M c=0or M d=0(see eq.[5]).Equation(6)implies that planets with comparable masses cannot have precessional degeneracy with g1=g2(Kinoshita&Nakai2000).The present value ofγd is0.29for planets c and d around Ups And.Today,in the absence of any residual disk(N c,d=0),g2∼g d<g1∼g c and the outer two planets of Ups And are not in a state of precessional degeneracy.Nevertheless,their present orbits can be approximated by equation(3)with the magnitude of A c,d much smaller than that of B c,d. Thus,the two planets primarily precess with eigenfrequency g1and their relative longitudes of periastrons librate over a restricted range of phases.This phase lock is equivalent of a state of secular resonance.This resonant interaction is the result of a dynamical feedback through the planets’secular interaction.2.4.Three Classes of Relative OrbitsIn order to illustrate the importance of this feedback effect,we substitute x≡e c/e d and η≡̟c−̟d,equations(1)and(2)reduce todx=(1−γd)+(C/x)(1−γd x2)cosη+∆N=A1+A2(x)cosη,(9) dτwhere∆N≡N c−N d,A1is given in equation(7),and A2(x)≡(C/x)(1−γd x2)is introduce for notational simplicity.There are three families of relative orbits which can be illustrated below with the linearized approximation solutions of equations(8)and(9)(see Appendix B).We show below that the relative magnitude of A1and A2determine the nature of the orbits.1)Circulation.In these solutions,the magnitude of x modulates about some values x0 such that x=x0+δx(τ)with an amplitude|δx(τ)|≪x0for all values ofηwhich ranges between0and2π.In the limit that|A1|≫|A2|,ηdecreases monotonically(because A1<0) whileδx oscillates(see solutions in eqs.[B7]and[B8]in Appendix B.2).To the lowest order these solutions reduce toCη≃A1τ,δx≃−Near the secular resonance where|A1|is relatively small,the amplitude ofδx becomes large.2)Libration.In the opposite limit that|A1|≪|A2|,the non axisymmetric secular interaction between the planets is important.There are stationary points in the(x-η)plane which center on the values of x=x m andη=0orπ(see Appendix B.1).Around these points,there are orbits with both small amplitude modulation such thatx−x mǫ≡x m 1+γd x2m (12) is the oscillation frequency(see Appendix B.1).The dimensionless amplitude of these libra-tional orbitsǫ0≪1.For these orbits,although g1=g2,secular interaction induces them to precess at similar frequencies that their relative longitude of periapse passage is always approximately aligned or anti aligned.Thus,these planets are effectively in a state of secular resonance.Note that because the libration is centered aroundη=0orπ,there is no effective angular momentum transfer despite the phase lock.3)Excitation.For systems with|A1|<|A2|,there are also orbits in which x andηhave very different values as x m and0(orπ)respectively.In these cases,ηwould evolve rapidly to a phase angleη1≃cos−1 −A1C(1−γd x2) ,(13) such that dη/dτis reduced to zero(see eq.[9])and the two planets become phase locked.For all non zero(orπ)values ofη1,the monotonic increases/decreases of x correspond to eccentricity excitation/damping,analogous to the situation of precessional degeneracy.In this state of near secular resonances,angular momentum is monotonically transferred from one planet to the other,resulting in a monotonic evolution of x.The modification of x in turn leads to an evolution inη.Along the path of x andη1evolution,g1=g2and the planets would librate about their evolving guiding center.For the special cases when the two planets enter the resonance,the second order solution reduces to that in equation(11)and all orbits become instantaneously librational,even thoughǫ0andη0may become arbitrarily large.At the center of resonance where A1,A2,and A1/A2all vanish,η1=π/2and x m evolves exponentially.We discuss the evolution of these systems in Appendix C.3.MODELSThe above analytic approximation is useful for isolating the three family of orbits along with two stationary points in the(x-η)diagram.These orbits generally follows the contours of equi Hamiltonian map(see Appendix D).However,the topological evolution of the Hamiltonian map alone is insufficient for the determination of the planets’orbits during the depletion of the disk.We carry out,below,numerical integration of the full equation of motion for the planets subject to the potential of their host star and nascent disks.In this section,we briefly describe a model prescription with which we examine the passage of librational degeneracy during the epoch of disk depletion.3.1.Planetary Formation Scenarios and Disk ModelWefirst discuss the physical process of disk depletion.According to conventional theo-ries,planets are formed through the condensation of grains which grow to planetesimals via cohesive collisions(Hayashi,Nakazawa,&Nakagawa1985;Lissauer1987;Wetherill1990). Upon attaining a sufficiently large mass,planetesimals accrete gas(Mizuno1980;Boden-heimer&Pollack1986;Pollack et al.1996).Eventually,their growth is terminated when protoplanets can tidally induce the formation of a gap near their orbit(Goldreich&Tremaine 1980;Lin&Papaloizou1980,1993;Takeuchi,Miyama,&Lin1996).Thereafter,gas in the inner of the disk continues to diffuse inward as it loses angular momentum to the planet.Since gas replenishment is cut offby the formation of the gap,the inner region is depleted well before the outer region of the disk.For the present discussions, we assume that the mass in the inner regions of the disk becomes negligible by the stage when the second planet is fully grown.Thus,in our analysis of planets’dynamical evolution, we only need to consider their interaction with the disk region extended outside the outer planet.Both the gravitational and tidal influences of the disk on the outer planet are much more intense than on the inner planet because the former is much closer to where most of the gas is distributed.When the planets’orbits lie in the same plane as their nascent disks and their distance to the nearest edge of their nascent disks is greater than the disks’scale height,a thin-disk approximation is sufficient for the computation of the disks’gravitational potential. Following the approach of Ward(1981),the self-gravitating potential for a disk model with a surface density profileΣ=Σ0(r0/r)k(withΣ0and r0being somefiducial values)can beexpressed asV(r)=2πGΣ0r r02n+k−1 rthe structure of the disk,we also consider a model in which the surface density of the diskis assumed to decline by an identical reduction factor everywhere such thatΣ0(t)=Σ0(t=0)exp(−t/t∆N)with r edge=constant.In the limit that a c/a d≪1and a c,d/r edge≪1,we find that the two prescriptions of the disk depletion do not lead to a significant difference inthe numerical results because∆N∝Σ/r5/2edge.For the approximate prescription of∆N(seeeq.[D2]),the magnitude of A1=1+∆N−γd changes on a time scaleτA1∼t∆N and the condition of the secular resonance A1=0is satisfied whenΣ0r3/204πM dn d −1.(15)This condition can be attained with both prescriptions.The orbital evolution of the planets can be directly integrated numerically,including the point-mass potential of the host stars and sibling planets as well as the contribution from the disk potential where the gravity can be expressed asF=dV4n+1A n r5.5AU.The panel(3)represents the epoch of secular resonance passage.At this stage,theprecession frequencies of the longitudes of periastrons nearly match with each other(i.e.,˙η∼0compared with˙x).Consequently,ηvaries slowly(withfinite values other than0orπ)in comparison with the non negligible evolution of x such that the eccentricities can changesignificantly.The equi-Hamiltonian map is locally symmetric about the lines ofη=−π,−π/2,0,π/2,andπand that of x=γ−1/2d∼2.When the nebula edge is beyond5.5AU,the precession speed of longitude of planet c’speriastron is faster than that of planet d’s periastron(except for the occasional regressions).Panel(4)shows the case that the nebula edge has retreated to7AU.The centers of theclosed libration tracks are reversed from the situation in panels(1)and(2).The planetaryconfiguration with e c>e d and̟c∼̟d in panel(2)becomes that with e c<e d and̟c∼̟d in panel(4).When the entire protoplanetary disk is depleted(panel(5)),stable orbits libratearound(x=0.5,η=0)or(x=7,η=±π).The shape of equi-Hamiltonian map is turned.Figure3is the time evolution of the equi-Hamiltonian map for upside down with x=γ−1/2dthe uniform-depletion prescription with r edge=4.5AU.In this model,the secular resonanceoccurs whenΣ0(t)/Σ0(t=0)=0.5.Although the time of the secular-resonance passageis slightly modified from that obtained with the model of inside-out depletion,the generalevolutionary pattern is not qualitatively changed.We also numerically integrate the orbits of planets c and d using the the same inside-outdisk depletion prescription as above(i.e.,we arbitrarily set the initial disk edge to be at4.2AU and specify its retreating speed to be10−5AU year−1).We consider two sets of initialconditions.Figures4a and c illustrate the evolution of a model with initial values of x=8andη=π/4.We illustrate the orbital evolution of planets c and d with the(x-η)diagram(see Fig.4a)and with e c and e d as a function of time(see Fig.4c).These results show thatthe planets’orbits initially librate relative to each other.The planets remain on libratingclosed tracks in the(x-η)diagram after the disk is totally depleted.But,as a consequenceof the secular-resonance passage,the eccentricity of planet c becomes smaller than that ofthe planet d.The time-averaged eccentricities after the disk depletion are e c ≃0.12and e d ≃0.2.We also consider a second set of initial conditions with x=8(e c=0.4,e d=0.05)andη=3π/2.In this case,the planets’orbits initially circulate relative to each other,but,theybecome trapped in librating orbits during the passage through the secular resonance(seeFig.4b).The mean eccentricities after the total depletion of the disk are very similar tothat obtained in the case of panel(a),because they are determined by the conservation ofangular momentum(eq.[D6]).We also carried out numerical integration of the full equations with the uniform-disk-depletion prescription with t∆N=105years,r edge=4.5AU.When the same initial conditionis used(i.e.x=8andη=π/4)were used,the pattern of the orbital evolution with thisprescription(see Fig.5a)is not significantly changed from that in Figure4a.In panelsb and c,we show the evolution from(x=8,η=π)and(x=5.6,η=0).The pair ofplanets in panel b has an initially wide open circulating track.During the disk depletion,they temporarily enter into closed librating track.But,their orbital configuration becomeswide open again when most of the disk material is depleted.The results of these numericalcalculations show that as long as the depletion time scale of the disk is longer than theoscillation period of eccentricity,the planets’orbits evolve adiabatically.Those systemswith librating orbits in(x-η)diagram prior to the disk depletion usually remain on the closedlibrating tracks after disk is totally depleted despite large changes in magnitude of x as in thecase of Ups And(see below).Those pairs of planets with the widely open circulating tracksinitially generally remain on open circulating tracks after the disk depletion is completedas in the case of HD168443(see below).Under some circumstances,it is also possible forplanets with marginally circulating/librating initial orbits to undergo transition after thedisk depletion.4.2.Orbital Evolution of Circulating Planets Around HD168443Next we consider the case of the planetary system around HD168443for which weassume the mass ratio to be M b/M c=0.451.Figure6shows the time evolution of the equi-Hamiltonian contour map of the HD168443system.The nebula edge retreats from insideto outside(from panel1to panel5).The secular resonance occurs when∆N=−0.86withr edge=7.2AU.Similar to the case of the planetary system around Upsilon Andromeda,thelibrating orbit with|η|<π/2are confined in regions with e b>e c prior to disk depletion.But,in the case of HD168443system,the domain of closed librating tracks is small.(Incontrast,for small M b/M c(a b/a c)1/2,the domain of open circulating tracks is large as isthe case for planets c and d around Ups And).The domain of closed librating tracks with|η|<π/2moves downward in the Hamiltonian map as disk depletion proceeds.If planets of HD168443initially follow open circulating tracks prior to the disk depletion,they wouldremain on the open circulating tracks after the disk depletion.In this case,during the onsetof secular resonance,planets b and c may briefly attain a librating track before moving ontoan open circulating track as the disk continues to deplete.Throughout the epoch of diskdepletion,the planets’longitudes of periastrons are widely separated such that little angularmomentum may be exchanged between them.Consequently,the net change of eccentricitiesis limited.。

a r X i v :a s t r o -p h /0702203v 1 7 F eb 20072Fathi Namounilarge amount of gas from the protoplanetary disk.Only after they formed,did they travel all the way to meet their current orbits.It is interesting to note that the concept of radial migration was already known in the contexts of accretion disks[32]in binary star systems and of planetary rings[18].Only before1995,one could not plausibly contemplate the prospect of suggesting the existence of massive planets that traveled all the way from Jupiter’s cur-rent location just to stop on a close orbit with a few-day period.The basic aspects of the process of planetary migration through the tidal interaction of a planet with the gaseous protoplanetary disk are now well understood[65]yet an important challenge remains:what stops planetary migration towards the star?The leading contender for stopping planetary migration is the planet’s interaction with the stellar magnetosphere but a definitive quantitative de-scription is still lacking.Extrasolar planetary eccentricities have equally resulted in a drastic change of perception:it is often heard that it is not the extrasolar planets that are eccentric,rather it is the solar system that lacks eccentricity.This perception is encouraged by the availability of some simple instabilities that one can set up in a many-body gravitational system to simulate the generation of the wild orbits of extrasolar planets.Upon close examination such instabilities as well as other eccentricity scenarios do not tell the whole story of how extrasolar planets become eccentric.In fact,just as the features of the planetary migra-tion process yield constraints on the planetary formation scenarios,so do the various theories of the eccentricity excitation.It is the aim of this chapter to review the various processes of the origin of extrasolar planets’eccentricities in the context of planetary formation.We start by reviewing the properties of extrasolar planetary orbits in section2. Section3contains a commentary on the various known theories of eccentricity excitation.Section4specializes in a recent addition to the eccentricity theories based on an relationship between the planets and the stellar jet that is powered by the protoplanetary accretion disk.Thefinal section5discusses how the eccentricity origin problem may contribute further to the theory of planet formation.2Eccentricity observationsExtrasolar planets are detected with various observational techniques[55]. The Doppler analysis of the reflex velocity of the host star is by far the most successful technique to date.It is also the technique that has uncovered the large eccentricities of extrasolar planets.If a planet has a circular orbit,the analyzed stellar spectrum yields a sinusoidal oscillation of the stellar reflex motion.If the planet is on an eccentric orbit,the reflex motion as a function of time becomes distorted with respect to a pure sine reflecting the unequal times the star spends in different locations along its orbit around the center of mass of the star-planet system(Figure1).The discovery of large eccentricityEccentric extrasolar planets3 orbits by the Doppler reflex velocity method is due to its ability to detect planets on wider orbits in contrast to that to the transit method.Planets on many-day periods have usually undergone tidal circularization by the host star.Fig.1.Reflex velocity of the stars HD75289(left)[63]and HD108147(right) [51].The planet around HD75289has a circular orbit while that around HD108147 has an orbital eccentricity of0.5.Pictures taken from the Geneva Extrasolar Planet Search http://obswww.unige.ch/˜u dry/planet/planet.htmlThe statistical analysis of extrasolar planet eccentricities reveals very few clues as to the origin of the elongated orbits.For the known sample of196 planets discovered by the reflex velocity and transit techniques,the median eccentricity is at0.21if all planets are counted and at0.28if planets with periods smaller than5days are excluded because their circular orbits simply reflect tidal circularization.The prevalence of such large eccentricities and the large typical mass of the detected planets(comparable to Jupiter’s)has encouraged the comparison of the extrasolar planetary systems to binary star systems.Depending on the methods used,similarities in the eccentricity dis-tribution of both populations can be found[57]or not[21].What is agreed upon is that there is no correlation between the size of the orbits and their eccentricities in each population,and no striking resemblance of the scatter of both populations in an orbital size versus eccentricity plane.The size of the orbit usually refers to either the semi-major axis or the pericentre radius.The latter is used to account for those orbits that have not yet had enough time to be circularized–as the pericentre distance is conserved under stellar tides. Finally,eccentricities show a vague correlation with the planetary masses with heavier planets enjoying larger eccentricities.4Fathi Namouni3Eccentricity origin theoriesSeven known explanations have been put forward to account for the large eccentricities of extrasolar planets.They are:(1)planet-planet scattering,(2) planet-protoplanetary disk interaction,(3)Kozai’s secular cycles,(4)excita-tion through radial migration into a mean motion resonance,(5)Stellar en-counters,(6)stellar-like N-body relaxation,and(7)excitation through stellar jet acceleration.In the following,we comment on these possibilities by dis-cussing their instability types,characteristic timescales,their epoch of appli-cability as well as their advantages and drawbacks.3.1Planet-planet scatteringPlanet-planet scattering is a simple process to generate eccentric orbits in an N-body gravitational system.If a system of two or more planets on planar cir-cular orbitsfind themselves“initially”closer than is permitted by Chirikov’s criterion for the overlapping of mean motion resonances[68],the planets scat-ter offone another leading to a system with more stable albeit eccentric config-urations.Depending on the number,masses and“initial”spacings of the plan-ets,the instability timescale varies between103to107years[52,67,15,40,16]. The epoch that is referred to by the adjective“initial”is that of the disap-pearance of the agent or the conditions that kept the planets from scattering offone another in thefirst place.This epoch is customarily associated with a significant dispersal of the parent gaseous protoplanetary disk.As well shall point out in the next section,planet-disk interaction is known to primarily erase orbital eccentricities.An additional condition for planet-planet scatter-ing to be operational is the absence of a significant population of smaller bodies such as the primordial asteroid belt.Depending on the mass spectrum in the planetary system,the smaller populations are able to limit the growth of the planetary eccentricities through dynamical friction[4].This at least how it is believed that the terrestrial planets in the solar system did not acquire large eccentricities[1,44].Numerical works that tackle the extrasolar eccen-tricity problem using planet-planet scattering do not consider the effect of leftover small-body populations after the gaseous disk has dispersed.Planets are set up at a few Hill radii from one another and initial conditions are sam-pled to reproduce the eccentricity of certain observed systems.The general excitation trend of planet-planet scattering leads to larger eccentricities than the ones observed.What may prove to be a serious problem for planet-planet scattering is the eccentricity distribution obtained in multiple systems that contain Jupiter-mass planets as well as Earth-mass planets.The conservation of angular momentum in this case will force the much smaller planets to have much larger eccentricities than the Jovian planets.Eccentric extrasolar planets53.2Planet-disk interactionA planet embedded in a gaseous disk excites sound waves at the locations of its mean motion resonances within the disk akin to the gravity waves excited by Saturn’s satellites in its ring system.The density enhancements at the mean motion resonances act back on the planet resulting in gravitational torque. Two types of resonances contribute to this torque:(1)corotation resonances that primarily affect the semi-major axis and tend to damp any acquired ec-centricity and(2)Lindblad resonances that primarily affect the eccentricity and tend to increase it[17,18].The torque contribution of the former is larger than the latter’s by about5%.Atfirst sight,planet-disk interaction damps the eccentricity on timescale that depends strongly on the disk’s thickness and less strongly on the disk’s mass density and the planet’s mass[64,2,37].The torques originating from higher order resonances as well as those pertaining to the relative inclination of the planet and the disk do not change the out-come significantly[19,45].Only if the corotation torque saturates,can the Lindblad resonances increase the eccentricity[19,20].The conditions under which saturation arises are difficult to quantify explaining why an eccentricity increase due to a disk-planet interaction has never been observed in numerical simulations although this might be due to numerical artefacts[38].3.3Secular Kozai cyclesIn his study of asteroids perturbed by Jupiter on high eccentricity and incli-nation orbits,Kozai[28]showed that the averaging of the interaction poten-tial over the mean motion without expanding the force amplitude in terms of eccentricity and inclination leads to new types of secular resonances.The conservation of the vertical component of angular momentum(vertical refers to the direction of Jupiter’s orbital normal)shows that when the orbital ec-centricity increases,the inclination decreases.In particular,if objects are set up on inclined but circular orbits,large eccentricities can be achieved as the inclination decreases in its motion around the secular resonance cycle.The application of the Kozai cycle to the eccentricities of extrasolar planets as-sumes that there is a binary star on a not-too-far inclined orbit that perturbs the planet that formed in a circular orbit in a timescale shorter than the Kozai libration cycle.In essence,the secular Kozai cycle idea transforms the eccentricity problem into an inclination problem.In this sense,the observed planets do not possess a proper eccentricity but one that is forced by the stel-lar binary and that will always oscillate between its original value,zero,and a maximum value depending on the planetary-binary semi-major axis ratio,the binary’s mass and its orbital inclination with respect to the plane on which the planet initially formed.When applied to specific binary star systems with one planet,the Kozai mechanism worksfine and helps characterize the orbit and mass of the secondary star required to excite eccentricity[27,42,14].Sta-tistically,Kozai based excitation of one-planet binary systems tend to yield6Fathi Namounilarger eccentricities than observed[58].As T Tauri stars form in multiple sys-tems,it is not unreasonable to try and apply the Kozai mechanism to the whole sample of observed extrasolar planets.The problem is that the Kozai cycle is usually destroyed by mutual gravitational interactions.The addition of more planets to the one-planet binary system,forces the precession of the planets’pericentres.If the planets are of comparable mass as it is observed in multiplanet systems,the Kozai cycle is lost.3.4Mean-motion resonancesThe role of mean motion resonances in exciting orbital eccentricity has its roots in the study of the orbital evolution of Jupiter’s and Saturn’s regular satellites under planetary and satellite tides[50].These satellite systems are known to be in or to have crossed mean motion resonances thereby acquiring forced eccentricities.The combined modeling of the orbital evolution,capture into resonance and the tidal interaction lends valuable bounds on the dissipa-tion factors of Jupiter,Saturn and their satellites.Extrasolar planets form in a gaseous disk that does not dissipate after they acquired most of their masses or else hot Jupiters would not exist.Planet-disk interaction naturally gives rise to orbital migration with different planets in the same system migrating at different rates.This differential migration makes planets in the same sys-tem encounter mean motion resonances.Capture into resonance may occur depending on whether the migration is convergent or divergent(for instance if the outer planet is moving faster or slower than the inner planet).Con-vergent migration leads to capture into resonance.The subsequent common migration of a planetary pair in resonance pumps up the eccentricities on the migration timescale[31,70].Divergent migration does not lead to resonance capture,instead eccentricity jumps are acquired at resonance passage[7,62]. While convergent migration is certainly the way the known resonant multiple systems have acquired their eccentricities,this excitation method involves a mystery that may shed light on how to halt planet migration in a gas disk. The mystery consists of the observation that convergent migration is far too efficient in exciting eccentricities to the point where in many systems,when capture occurs,migration must stop quickly thereafter or else eccentricities are pumped up to much larger values than those observed.As it is implausi-ble to invoke the dispersal of the gas disk,planetary migration may become ineffective because of the nonlinear response of the gas disk to the planet pair.It is interesting to note that when capture occurs,the planetary relative inclinations may be excited as unlike planetary satellites that orbit Jupiter and Saturn,the central potential is keplerian.Consequently,for the same or-der,eccentricity and inclination resonances are close(but not coincident as the gas disk modifies the pericentre and node precession rates).Planet-disk interaction is not well understood for large eccentricity planets and off-plane (inclined with respect to the disk)orbits.The often used formulas for eccen-tricity damping from the disk torques have not been verified for eccentric andEccentric extrasolar planets7 inclined planets.Divergent migration has the advantage of being applicable to the wider non-resonant multiplanet systems.For a planetary pair,divergent migration requires that the inner planet migrates faster than the outer one. Gap driven migration(also known as type II)is favorable to such a condition as the migration rate is the viscous timescale of the disk.Divergent migration may take place because viscosity is likely to be a decreasing function of the distance to the star.If however the part of the disk that is located between the two planets is dispersed as when the two planetary gaps merge,the direction and rate of migration may be altered significantly[39].3.5Stellar encountersStellar encounters are common events in star clusters.A planetary system that encounters a star will have its planets feel a tidal force that elongates their orbits.For inner planets that orbit close to the host star,the excitation which lasts for about1000years will occur on a secular timescale.Outer planets if they exist will feel a localized impulse somewhere in their orbits.Typical encounter frequencies are of one in5×109years while typical encounter parameters are a few hundred AU.Unless planets are way outside the classical planetary region(inside30AU),excitation is not efficient[71].To reverse this conclusion and account for the eccentricities of inner planets,the system must contain several planets with increasing distance and mass from the star in order propagate the stellar tug felt by the outermost planet down to the innermost ones[71].3.6Stellar-like relaxationThe qualitative similarity of the eccentricities of extrasolar planets and stellar binaries suggests that planets may form through similar processes as those of multiple stellar systems.If planets formed by gravitational instability,the for-mation time is so short that the planetsfind themselves confined to a smaller space than their orbital stability permits.The relaxation of such systems leads to some planetary ejections and many large eccentricity orbits[49].The ap-plicability of this scenario is limited by two facts:first,the minimum planet mass the gravitational instability allows is a few Jupiter masses.This means that stellar-like relaxation does not work for planets with masses comparable to or smaller than Jupiter’s.Second,if a two-phase formation where small planets form through rocky core accretion and the larger ones through gravi-tational instability[60],then it is likely that the relaxation of the larger plan-ets destroys the smaller planets.This is because the gravitational instability timescales are usually smaller than the planetesimal accumulation timescales. In fact,if large mass planets form through gravitational instability,they are likely to inhibit planetesimal accumulation by clearing the inner disk before planetary embryos are born.8Fathi Namouni4Jet-induced excitationStellar jets enter the eccentricity excitation problem because of their ubiquity and simplicity[46].The story of how this works is as follows:although there is disagreement on whether there is a statistically significant resemblance be-tween the eccentricity distributions of extrasolar planets and stellar binary systems,the qualitative similarity is beyond doubt.Those who wish for the similarity to be quantitative,would like to affirm the view that planets are the lower end of the outcome of star formation.This question has been settled observationally in2005with two observations:thefirst is a hot Saturn with a giant rocky core discovered by combining the Doppler reflex velocity method with transit photometry[54].The second is the imaging of thefirst planetary candidate which because of bias due to contrast and resolution happens to be a warm distant companion orbiting a young brown dwarf[5].This proves that planets do not need large rocky cores and may form by gravitational instability.Exit the link between how planets form and their eccentricities.If planets do not form like binary stars,perhaps they undergo similar ex-citation processes that lend them similarly elongated orbits.In view of the different physical environments where planets and stars form,the simplest possible excitation process may depend weakly or not at all on the local dy-namics of the stellar or planetary companion.Mathematically,this amounts to saying that the acceleration imparted by the process is independent of position and velocity.Simplicity therefore dictates that the process imparts a constant accelera-tion that operates during afinite time window.Simplicity also comes with two added advantages:we can already know the excitation time scale and the min-imal acceleration amplitude.Dimensional analysis shows that the excitation timescale has to be proportional to v/A where v is the keplerian velocity of the companion around the main star.Further,if the acceleration is to achieve its purpose within the lifetime of the system,v/A must be smaller than about109 years.This tells us that the acceleration A>3×10−16(v/10km s−1)km s−2.The process lacks one more attribute:direction.If the acceleration is inde-pendent of the formation processes,its direction cannot depend on anything related to the planetary companion such as its orbital plane or the direction from the star to the companion.In an inertial frame related to the planetary or stellar system,we are not left with much choice but the star’s rotation axis.To sum up,what we are looking for is a process that appears everywhere where planet and star formation takes place,acts like a rocket(i.e.with an acceleration that does not depend on the position and velocity of the system) and whose direction is related to the star’s rotation axis.The answer is stellar jets[12,23].Do planets exist when jets are active?The answer is quite likely.Known hot Jupiters have moved close to their host stars because of their interaction with the gas disk.So we know the gas disk was present and had viscosity well after planetsfinished forming.The gas accreting on the star because of viscosityEccentric extrasolar planets9 is the main ingredient along with the magneticfield that threads it needed to launch stellar jets and disk winds.It would therefore be an interesting coincidence that jets shut offwhen planets appear in the gas disk a few AU away from the star well outside the jet launching region.Are there any observational hints that jet-sustaining disks contain planets? The only possible hint so far is the observation of variable brightness asym-metries in some jet-sustaining disks[8,56,66].The variability timescales of a few days to a few years are so small that they imply either a peculiar stellar activity in the form of single hot spots or the presence of distortions in the disk at the location where the orbital period matches the variablity timescale. Thefirst option requires a complex stellar magneticfield that differentiates strongly between the two stellar poles.The second option may be caused dis-tortions in the disk whose origin could be the presence of embedded compact objects.Do jets have enough strength to build eccentricity?Inferred mass loss rates for known young T Tauri stars lie in the range∼10−8M⊙year−1to 10−10M⊙year−1and may be two orders of magnitude larger depending on the way the rate is measured from the luminosity of forbidden lines[22,25,29]. The jet also needs to be asymmetric with respect to the star’s equator plane or else there would be no acceleration.Interestingly,a growing number of bipolar jets from young stars[24,13,30,35]are known to be asymmetric as the velocities of the jet and counterjet differ by about a factor of2.Mass loss processes in young stars therefore yield accelerations:A∼10−13 ˙M300km s−1 M⊙10Fathi Namouniwhere v g is the velocity dispersion in the Galaxy yields:τ≤1053×10−12km s−210km s−1years.(2)This timescale is shorter than the disk’s lifetime.In practice,we shall see that shorter times are needed.Further excitation properties can be deduced by analyzing the effect of the combined jet-induced acceleration and the star’s gravitational attraction.As the star’s pull decreases with distance,there is a specific location where the latter matches the jet-induced acceleration(that is independent of position and velocity).Outside this radius,the star’s pull is weak and orbits escape its gravity.This reveals an interesting feature of jet-induced acceleration:stellar jets are responsible for the outer truncation of circumstellar disks.It is clear that in the interior vicinity of the truncation radius,the orbital perturbations are large as the excitation time becomes comparable to the orbital period.In this region,the keplerian orbits are subject to a sudden excitation;not only the eccentricities are excited but the semi-major axes are also affected leading to inward or outward migration.Well inside the truncation radius,the excitation time is much larger than the orbital period.In this region,eccentricity builds up slowly over a large number of revolutions of the planet around the star and the mean orbital radius remains constant on average.Excitation in this region occurs on secular timescales.Planetary companions mostly fall inside the secular region as they are far inside the truncation radius which is more or less the size of the protoplanetary disk.4.1Secular jet-induced excitationIn the secular region where the excitation time is larger than the orbital period, the dynamics of excitation can be simplified by averaging the acceleration over the orbital period of the companion.For a constant acceleration,the interaction potential is simply R=A·x where x is the position vector. Averaging the interaction potential amounts to averaging the position vector of a pure keplerian motion.A simple calculation shows that x =−3ae x(f= 0)/2r where f,a,e and r are the true anomaly,the semi-major axis,the eccentricity and radius of the keplerian orbit.The direction of x at pericentre is that of the eccentricity vector e=v×h/G(m+M)−x/|x|where v is the velocity vector of the companion,G is the gravitational constant,m and M are the masses of the companion and the host star and h=x×v is the specific angular momentum.This enables us to write the secular potential as: R =−32A ae sin(̟−Ω)sin I,(3) where in the last equality,the z–direction of the reference frame is chosen along along A and̟,Ω,I,are the longitude of pericentre,longitude of ascendingnode and the inclination of the orbit.To simplify the excitation problem further,we use the conservation of the component of angular momentum h along the direction of acceleration as A ·˙h=A ·(x ×A )=0.In the reference frame where A is along the z –direction,the conservation of angular momentum yields (1−e 2)1/2cos I =cos I 0where I 0is the initial inclination of the keplerian orbit with respect to the jet-induced acceleration.This relation enables us to eliminate the inclination variable in R and reduce the problem to an integrable,one dimensional system with:R =−3sin 2I 0−e 21−e 2∂ω,˙ω=√na 2e ∂ RG (M +m )/a 3is the companion’s mean motion.In this one-dimensional system,e and ωfollow curves of constant R shown in Figure (2).There are equilibria at ω=±90◦and e =√cos I 0).The maximum value of e is sin I 0and correspondsto the cycle of initially circular orbits.For these orbits, R =0throughout their cycle implying that the orbits orientation can take only one value ω=0modulo 180◦.-180-90090180w 00.10.20.30.40.5eFig.2.Contour plots of the acceleration potential (4)in the eccentricity e and argument of pericentre ω(◦)plane.The direction of acceleration makes an angle I 0=30◦with respect to the companion’s angular momentum vector.The time evolution of the two orbits (e =0,ω=0)and (e =0.3,ω=90◦)is shown in Figure(3).For time-dependent accelerations and provided that the variation timescale is longer than the orbital period,the eccentricity evolution is given by:˙e=3A(t)ǫsin2I0−e2,(6)whereǫis the sign of cosωwhich is set by the requirement that e≥0.The solution of(6)can be found exactly as:e(T)= sin 31−e(T)2.For strictly con-stant accelerations(infinite time window),A(t)=A0and e oscillates between 0and sin I0at the excitation frequency:n A=3|A0|0 0.10.20.30.40.5 10000 20000 30000 40000 50000t e 09018010000 20000 30000 40000 50000t w-30 030 10000 20000 30000 40000 50000t IFig.3.Time evolution under a conservative acceleration.The eccentricity e ,ar-gument of pericentre ω(◦)and inclination I (◦)are shown for an initially circular orbit e =0(solid)and an orbit librating about the secular resonance ω=90◦with an initial eccentricity e =0.3(dashed).The semi-major axis is identical for both orbits and is set to unity.The acceleration corresponds to a period of 104years at 1AU.The plots were obtained by the numerical integration of the full equations of motion.disk’s plane resulting from interactions with stellar companions as T Tauri stars are known to form in multiple systems.Precession is attractive because it offers the possibility of resonance if the excitation frequency n A matches the jet precession frequency ΩA .This in fact is exactly what happens when the eccentricity evolution is derived in the situation where the constant magnitude acceleration rotates at a constant rate.It turns out that the corresponding secular problem is also integrable.The eccentricity and inclination evolution are given by [46]:e 2=p 2sin 2α0 0.10.20.30.4 0.50 10000 20000 30000 40000 50000t e Fig.4.Eccentricity excitation by time-dependent constant-direction accelerations.The equations of motion are integrated numerically with an acceleration A (t )=A 0H (t )exp −t/τwhere H (t )is the Heaviside unit step function and A 0=2.21×10−11km s −2.The oscillation period at 1AU is 1.11×105years.The timescale τ=7200is chosen so that V =5km s −1.The curves correspond the semi-major axes:1AU (solid),32AU (dashed)and 128AU (dotted).cos I =11−e 2(p 4−p 2+2+p 2[p 2−3]cos 2α)+p 2sin α2(p 2+1+2p cos α)cos ν+t (10)+p 2sin α2(p 2+1−2p cos α)cos ν−t (11)where αis the jet angle with respect to the z –axis of the reference frame,ν2±=p 2+1∓2p cos α,p =n A /2ΩA ,and the time t is normalized by ΩA .The companion’s initial orbit is circular and lies in the xy –plane.Nominal resonance is defined where the frequency match,p =1,occurs.It corresponds to a nominal resonant semi-major axis a res given as:a res ≃4 M +m2×10−10km s −2 −2 T prec a/a res .Far inside resonance (p ≪1),the jetprecesses faster than the eccentricity excitation leading to a reduction of the eccentricity amplitude from sin αto 2p sin α.Far outside resonance (p ≫1),the jet’s precession is slow compared to the eccentricity excitation so that the latter is described by a constant acceleration without rotation.In the resonance region,the proximity of p to unity increases the denominators of the eccentricity expression (9)which leads to eccentricities close to unity.At exact resonance,the eccentricity reaches unity regardless of the jet angle.The width of the region around resonance increases with the jet angle α.These features are illustrated in Figures (5)where we plot the expressions (9)and(11)for a jet angle α=1◦,an excitation time 2π/n A =104years,and the three values of p :0.05,0.9,and 1.Finally,we note that as the eccentricity excitation time is n A ,no resonant forcing occurs when ΩA =n in the secular region (n A ≪n ).。

Observation of Gravitational Waves from a Binary Black Hole MergerB.P.Abbott et al.*(LIGO Scientific Collaboration and Virgo Collaboration)(Received21January2016;published11February2016)On September14,2015at09:50:45UTC the two detectors of the Laser Interferometer Gravitational-Wave Observatory simultaneously observed a transient gravitational-wave signal.The signal sweeps upwards in frequency from35to250Hz with a peak gravitational-wave strain of1.0×10−21.It matches the waveform predicted by general relativity for the inspiral and merger of a pair of black holes and the ringdown of the resulting single black hole.The signal was observed with a matched-filter signal-to-noise ratio of24and a false alarm rate estimated to be less than1event per203000years,equivalent to a significance greaterthan5.1σ.The source lies at a luminosity distance of410þ160−180Mpc corresponding to a redshift z¼0.09þ0.03−0.04.In the source frame,the initial black hole masses are36þ5−4M⊙and29þ4−4M⊙,and the final black hole mass is62þ4−4M⊙,with3.0þ0.5−0.5M⊙c2radiated in gravitational waves.All uncertainties define90%credible intervals.These observations demonstrate the existence of binary stellar-mass black hole systems.This is the first direct detection of gravitational waves and the first observation of a binary black hole merger.DOI:10.1103/PhysRevLett.116.061102I.INTRODUCTIONIn1916,the year after the final formulation of the field equations of general relativity,Albert Einstein predicted the existence of gravitational waves.He found that the linearized weak-field equations had wave solutions: transverse waves of spatial strain that travel at the speed of light,generated by time variations of the mass quadrupole moment of the source[1,2].Einstein understood that gravitational-wave amplitudes would be remarkably small;moreover,until the Chapel Hill conference in 1957there was significant debate about the physical reality of gravitational waves[3].Also in1916,Schwarzschild published a solution for the field equations[4]that was later understood to describe a black hole[5,6],and in1963Kerr generalized the solution to rotating black holes[7].Starting in the1970s theoretical work led to the understanding of black hole quasinormal modes[8–10],and in the1990s higher-order post-Newtonian calculations[11]preceded extensive analytical studies of relativistic two-body dynamics[12,13].These advances,together with numerical relativity breakthroughs in the past decade[14–16],have enabled modeling of binary black hole mergers and accurate predictions of their gravitational waveforms.While numerous black hole candidates have now been identified through electromag-netic observations[17–19],black hole mergers have not previously been observed.The discovery of the binary pulsar system PSR B1913þ16 by Hulse and Taylor[20]and subsequent observations of its energy loss by Taylor and Weisberg[21]demonstrated the existence of gravitational waves.This discovery, along with emerging astrophysical understanding[22], led to the recognition that direct observations of the amplitude and phase of gravitational waves would enable studies of additional relativistic systems and provide new tests of general relativity,especially in the dynamic strong-field regime.Experiments to detect gravitational waves began with Weber and his resonant mass detectors in the1960s[23], followed by an international network of cryogenic reso-nant detectors[24].Interferometric detectors were first suggested in the early1960s[25]and the1970s[26].A study of the noise and performance of such detectors[27], and further concepts to improve them[28],led to proposals for long-baseline broadband laser interferome-ters with the potential for significantly increased sensi-tivity[29–32].By the early2000s,a set of initial detectors was completed,including TAMA300in Japan,GEO600 in Germany,the Laser Interferometer Gravitational-Wave Observatory(LIGO)in the United States,and Virgo in binations of these detectors made joint obser-vations from2002through2011,setting upper limits on a variety of gravitational-wave sources while evolving into a global network.In2015,Advanced LIGO became the first of a significantly more sensitive network of advanced detectors to begin observations[33–36].A century after the fundamental predictions of Einstein and Schwarzschild,we report the first direct detection of gravitational waves and the first direct observation of a binary black hole system merging to form a single black hole.Our observations provide unique access to the*Full author list given at the end of the article.Published by the American Physical Society under the terms of the Creative Commons Attribution3.0License.Further distri-bution of this work must maintain attribution to the author(s)and the published article’s title,journal citation,and DOI.properties of space-time in the strong-field,high-velocity regime and confirm predictions of general relativity for the nonlinear dynamics of highly disturbed black holes.II.OBSERVATIONOn September14,2015at09:50:45UTC,the LIGO Hanford,W A,and Livingston,LA,observatories detected the coincident signal GW150914shown in Fig.1.The initial detection was made by low-latency searches for generic gravitational-wave transients[41]and was reported within three minutes of data acquisition[43].Subsequently, matched-filter analyses that use relativistic models of com-pact binary waveforms[44]recovered GW150914as the most significant event from each detector for the observa-tions reported here.Occurring within the10-msintersite FIG.1.The gravitational-wave event GW150914observed by the LIGO Hanford(H1,left column panels)and Livingston(L1,rightcolumn panels)detectors.Times are shown relative to September14,2015at09:50:45UTC.For visualization,all time series are filtered with a35–350Hz bandpass filter to suppress large fluctuations outside the detectors’most sensitive frequency band,and band-reject filters to remove the strong instrumental spectral lines seen in the Fig.3spectra.Top row,left:H1strain.Top row,right:L1strain.GW150914arrived first at L1and6.9þ0.5−0.4ms later at H1;for a visual comparison,the H1data are also shown,shifted in time by this amount and inverted(to account for the detectors’relative orientations).Second row:Gravitational-wave strain projected onto each detector in the35–350Hz band.Solid lines show a numerical relativity waveform for a system with parameters consistent with those recovered from GW150914[37,38]confirmed to99.9%by an independent calculation based on[15].Shaded areas show90%credible regions for two independent waveform reconstructions.One(dark gray)models the signal using binary black hole template waveforms [39].The other(light gray)does not use an astrophysical model,but instead calculates the strain signal as a linear combination of sine-Gaussian wavelets[40,41].These reconstructions have a94%overlap,as shown in[39].Third row:Residuals after subtracting the filtered numerical relativity waveform from the filtered detector time series.Bottom row:A time-frequency representation[42]of the strain data,showing the signal frequency increasing over time.propagation time,the events have a combined signal-to-noise ratio(SNR)of24[45].Only the LIGO detectors were observing at the time of GW150914.The Virgo detector was being upgraded, and GEO600,though not sufficiently sensitive to detect this event,was operating but not in observational mode.With only two detectors the source position is primarily determined by the relative arrival time and localized to an area of approximately600deg2(90% credible region)[39,46].The basic features of GW150914point to it being produced by the coalescence of two black holes—i.e., their orbital inspiral and merger,and subsequent final black hole ringdown.Over0.2s,the signal increases in frequency and amplitude in about8cycles from35to150Hz,where the amplitude reaches a maximum.The most plausible explanation for this evolution is the inspiral of two orbiting masses,m1and m2,due to gravitational-wave emission.At the lower frequencies,such evolution is characterized by the chirp mass[11]M¼ðm1m2Þ3=5121=5¼c3G596π−8=3f−11=3_f3=5;where f and_f are the observed frequency and its time derivative and G and c are the gravitational constant and speed of light.Estimating f and_f from the data in Fig.1, we obtain a chirp mass of M≃30M⊙,implying that the total mass M¼m1þm2is≳70M⊙in the detector frame. This bounds the sum of the Schwarzschild radii of thebinary components to2GM=c2≳210km.To reach an orbital frequency of75Hz(half the gravitational-wave frequency)the objects must have been very close and very compact;equal Newtonian point masses orbiting at this frequency would be only≃350km apart.A pair of neutron stars,while compact,would not have the required mass,while a black hole neutron star binary with the deduced chirp mass would have a very large total mass, and would thus merge at much lower frequency.This leaves black holes as the only known objects compact enough to reach an orbital frequency of75Hz without contact.Furthermore,the decay of the waveform after it peaks is consistent with the damped oscillations of a black hole relaxing to a final stationary Kerr configuration. Below,we present a general-relativistic analysis of GW150914;Fig.2shows the calculated waveform using the resulting source parameters.III.DETECTORSGravitational-wave astronomy exploits multiple,widely separated detectors to distinguish gravitational waves from local instrumental and environmental noise,to provide source sky localization,and to measure wave polarizations. The LIGO sites each operate a single Advanced LIGO detector[33],a modified Michelson interferometer(see Fig.3)that measures gravitational-wave strain as a differ-ence in length of its orthogonal arms.Each arm is formed by two mirrors,acting as test masses,separated by L x¼L y¼L¼4km.A passing gravitational wave effec-tively alters the arm lengths such that the measured difference isΔLðtÞ¼δL x−δL y¼hðtÞL,where h is the gravitational-wave strain amplitude projected onto the detector.This differential length variation alters the phase difference between the two light fields returning to the beam splitter,transmitting an optical signal proportional to the gravitational-wave strain to the output photodetector. To achieve sufficient sensitivity to measure gravitational waves,the detectors include several enhancements to the basic Michelson interferometer.First,each arm contains a resonant optical cavity,formed by its two test mass mirrors, that multiplies the effect of a gravitational wave on the light phase by a factor of300[48].Second,a partially trans-missive power-recycling mirror at the input provides addi-tional resonant buildup of the laser light in the interferometer as a whole[49,50]:20W of laser input is increased to700W incident on the beam splitter,which is further increased to 100kW circulating in each arm cavity.Third,a partially transmissive signal-recycling mirror at the outputoptimizes FIG. 2.Top:Estimated gravitational-wave strain amplitude from GW150914projected onto H1.This shows the full bandwidth of the waveforms,without the filtering used for Fig.1. The inset images show numerical relativity models of the black hole horizons as the black holes coalesce.Bottom:The Keplerian effective black hole separation in units of Schwarzschild radii (R S¼2GM=c2)and the effective relative velocity given by the post-Newtonian parameter v=c¼ðGMπf=c3Þ1=3,where f is the gravitational-wave frequency calculated with numerical relativity and M is the total mass(value from Table I).the gravitational-wave signal extraction by broadening the bandwidth of the arm cavities [51,52].The interferometer is illuminated with a 1064-nm wavelength Nd:Y AG laser,stabilized in amplitude,frequency,and beam geometry [53,54].The gravitational-wave signal is extracted at the output port using a homodyne readout [55].These interferometry techniques are designed to maxi-mize the conversion of strain to optical signal,thereby minimizing the impact of photon shot noise (the principal noise at high frequencies).High strain sensitivity also requires that the test masses have low displacement noise,which is achieved by isolating them from seismic noise (low frequencies)and designing them to have low thermal noise (intermediate frequencies).Each test mass is suspended as the final stage of a quadruple-pendulum system [56],supported by an active seismic isolation platform [57].These systems collectively provide more than 10orders of magnitude of isolation from ground motion for frequen-cies above 10Hz.Thermal noise is minimized by using low-mechanical-loss materials in the test masses and their suspensions:the test masses are 40-kg fused silica substrates with low-loss dielectric optical coatings [58,59],and are suspended with fused silica fibers from the stage above [60].To minimize additional noise sources,all components other than the laser source are mounted on vibration isolation stages in ultrahigh vacuum.To reduce optical phase fluctuations caused by Rayleigh scattering,the pressure in the 1.2-m diameter tubes containing the arm-cavity beams is maintained below 1μPa.Servo controls are used to hold the arm cavities on resonance [61]and maintain proper alignment of the optical components [62].The detector output is calibrated in strain by measuring its response to test mass motion induced by photon pressure from a modulated calibration laser beam [63].The calibration is established to an uncertainty (1σ)of less than 10%in amplitude and 10degrees in phase,and is continuously monitored with calibration laser excitations at selected frequencies.Two alternative methods are used to validate the absolute calibration,one referenced to the main laser wavelength and the other to a radio-frequencyoscillator(a)FIG.3.Simplified diagram of an Advanced LIGO detector (not to scale).A gravitational wave propagating orthogonally to the detector plane and linearly polarized parallel to the 4-km optical cavities will have the effect of lengthening one 4-km arm and shortening the other during one half-cycle of the wave;these length changes are reversed during the other half-cycle.The output photodetector records these differential cavity length variations.While a detector ’s directional response is maximal for this case,it is still significant for most other angles of incidence or polarizations (gravitational waves propagate freely through the Earth).Inset (a):Location and orientation of the LIGO detectors at Hanford,WA (H1)and Livingston,LA (L1).Inset (b):The instrument noise for each detector near the time of the signal detection;this is an amplitude spectral density,expressed in terms of equivalent gravitational-wave strain amplitude.The sensitivity is limited by photon shot noise at frequencies above 150Hz,and by a superposition of other noise sources at lower frequencies [47].Narrow-band features include calibration lines (33–38,330,and 1080Hz),vibrational modes of suspension fibers (500Hz and harmonics),and 60Hz electric power grid harmonics.[64].Additionally,the detector response to gravitational waves is tested by injecting simulated waveforms with the calibration laser.To monitor environmental disturbances and their influ-ence on the detectors,each observatory site is equipped with an array of sensors:seismometers,accelerometers, microphones,magnetometers,radio receivers,weather sensors,ac-power line monitors,and a cosmic-ray detector [65].Another∼105channels record the interferometer’s operating point and the state of the control systems.Data collection is synchronized to Global Positioning System (GPS)time to better than10μs[66].Timing accuracy is verified with an atomic clock and a secondary GPS receiver at each observatory site.In their most sensitive band,100–300Hz,the current LIGO detectors are3to5times more sensitive to strain than initial LIGO[67];at lower frequencies,the improvement is even greater,with more than ten times better sensitivity below60Hz.Because the detectors respond proportionally to gravitational-wave amplitude,at low redshift the volume of space to which they are sensitive increases as the cube of strain sensitivity.For binary black holes with masses similar to GW150914,the space-time volume surveyed by the observations reported here surpasses previous obser-vations by an order of magnitude[68].IV.DETECTOR VALIDATIONBoth detectors were in steady state operation for several hours around GW150914.All performance measures,in particular their average sensitivity and transient noise behavior,were typical of the full analysis period[69,70]. Exhaustive investigations of instrumental and environ-mental disturbances were performed,giving no evidence to suggest that GW150914could be an instrumental artifact [69].The detectors’susceptibility to environmental disturb-ances was quantified by measuring their response to spe-cially generated magnetic,radio-frequency,acoustic,and vibration excitations.These tests indicated that any external disturbance large enough to have caused the observed signal would have been clearly recorded by the array of environ-mental sensors.None of the environmental sensors recorded any disturbances that evolved in time and frequency like GW150914,and all environmental fluctuations during the second that contained GW150914were too small to account for more than6%of its strain amplitude.Special care was taken to search for long-range correlated disturbances that might produce nearly simultaneous signals at the two sites. No significant disturbances were found.The detector strain data exhibit non-Gaussian noise transients that arise from a variety of instrumental mecha-nisms.Many have distinct signatures,visible in auxiliary data channels that are not sensitive to gravitational waves; such instrumental transients are removed from our analyses [69].Any instrumental transients that remain in the data are accounted for in the estimated detector backgrounds described below.There is no evidence for instrumental transients that are temporally correlated between the two detectors.V.SEARCHESWe present the analysis of16days of coincident observations between the two LIGO detectors from September12to October20,2015.This is a subset of the data from Advanced LIGO’s first observational period that ended on January12,2016.GW150914is confidently detected by two different types of searches.One aims to recover signals from the coalescence of compact objects,using optimal matched filtering with waveforms predicted by general relativity. The other search targets a broad range of generic transient signals,with minimal assumptions about waveforms.These searches use independent methods,and their response to detector noise consists of different,uncorrelated,events. However,strong signals from binary black hole mergers are expected to be detected by both searches.Each search identifies candidate events that are detected at both observatories consistent with the intersite propa-gation time.Events are assigned a detection-statistic value that ranks their likelihood of being a gravitational-wave signal.The significance of a candidate event is determined by the search background—the rate at which detector noise produces events with a detection-statistic value equal to or higher than the candidate event.Estimating this back-ground is challenging for two reasons:the detector noise is nonstationary and non-Gaussian,so its properties must be empirically determined;and it is not possible to shield the detector from gravitational waves to directly measure a signal-free background.The specific procedure used to estimate the background is slightly different for the two searches,but both use a time-shift technique:the time stamps of one detector’s data are artificially shifted by an offset that is large compared to the intersite propagation time,and a new set of events is produced based on this time-shifted data set.For instrumental noise that is uncor-related between detectors this is an effective way to estimate the background.In this process a gravitational-wave signal in one detector may coincide with time-shifted noise transients in the other detector,thereby contributing to the background estimate.This leads to an overestimate of the noise background and therefore to a more conservative assessment of the significance of candidate events.The characteristics of non-Gaussian noise vary between different time-frequency regions.This means that the search backgrounds are not uniform across the space of signals being searched.To maximize sensitivity and provide a better estimate of event significance,the searches sort both their background estimates and their event candidates into differ-ent classes according to their time-frequency morphology. The significance of a candidate event is measured against the background of its class.To account for having searchedmultiple classes,this significance is decreased by a trials factor equal to the number of classes [71].A.Generic transient searchDesigned to operate without a specific waveform model,this search identifies coincident excess power in time-frequency representations of the detector strain data [43,72],for signal frequencies up to 1kHz and durations up to a few seconds.The search reconstructs signal waveforms consistent with a common gravitational-wave signal in both detectors using a multidetector maximum likelihood method.Each event is ranked according to the detection statistic ηc ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi2E c =ð1þE n =E c Þp ,where E c is the dimensionless coherent signal energy obtained by cross-correlating the two reconstructed waveforms,and E n is the dimensionless residual noise energy after the reconstructed signal is subtracted from the data.The statistic ηc thus quantifies the SNR of the event and the consistency of the data between the two detectors.Based on their time-frequency morphology,the events are divided into three mutually exclusive search classes,as described in [41]:events with time-frequency morphology of known populations of noise transients (class C1),events with frequency that increases with time (class C3),and all remaining events (class C2).Detected with ηc ¼20.0,GW150914is the strongest event of the entire search.Consistent with its coalescence signal signature,it is found in the search class C3of events with increasing time-frequency evolution.Measured on a background equivalent to over 67400years of data and including a trials factor of 3to account for the search classes,its false alarm rate is lower than 1in 22500years.This corresponds to a probability <2×10−6of observing one or more noise events as strong as GW150914during the analysis time,equivalent to 4.6σ.The left panel of Fig.4shows the C3class results and background.The selection criteria that define the search class C3reduce the background by introducing a constraint on the signal morphology.In order to illustrate the significance of GW150914against a background of events with arbitrary shapes,we also show the results of a search that uses the same set of events as the one described above but without this constraint.Specifically,we use only two search classes:the C1class and the union of C2and C3classes (C 2þC 3).In this two-class search the GW150914event is found in the C 2þC 3class.The left panel of Fig.4shows the C 2þC 3class results and background.In the background of this class there are four events with ηc ≥32.1,yielding a false alarm rate for GW150914of 1in 8400years.This corresponds to a false alarm probability of 5×10−6equivalent to 4.4σ.FIG.4.Search results from the generic transient search (left)and the binary coalescence search (right).These histograms show the number of candidate events (orange markers)and the mean number of background events (black lines)in the search class where GW150914was found as a function of the search detection statistic and with a bin width of 0.2.The scales on the top give the significance of an event in Gaussian standard deviations based on the corresponding noise background.The significance of GW150914is greater than 5.1σand 4.6σfor the binary coalescence and the generic transient searches,respectively.Left:Along with the primary search (C3)we also show the results (blue markers)and background (green curve)for an alternative search that treats events independently of their frequency evolution (C 2þC 3).The classes C2and C3are defined in the text.Right:The tail in the black-line background of the binary coalescence search is due to random coincidences of GW150914in one detector with noise in the other detector.(This type of event is practically absent in the generic transient search background because they do not pass the time-frequency consistency requirements used in that search.)The purple curve is the background excluding those coincidences,which is used to assess the significance of the second strongest event.For robustness and validation,we also use other generic transient search algorithms[41].A different search[73]and a parameter estimation follow-up[74]detected GW150914 with consistent significance and signal parameters.B.Binary coalescence searchThis search targets gravitational-wave emission from binary systems with individual masses from1to99M⊙, total mass less than100M⊙,and dimensionless spins up to 0.99[44].To model systems with total mass larger than 4M⊙,we use the effective-one-body formalism[75],whichcombines results from the post-Newtonian approach [11,76]with results from black hole perturbation theory and numerical relativity.The waveform model[77,78] assumes that the spins of the merging objects are alignedwith the orbital angular momentum,but the resultingtemplates can,nonetheless,effectively recover systemswith misaligned spins in the parameter region ofGW150914[44].Approximately250000template wave-forms are used to cover this parameter space.The search calculates the matched-filter signal-to-noiseratioρðtÞfor each template in each detector and identifiesmaxima ofρðtÞwith respect to the time of arrival of the signal[79–81].For each maximum we calculate a chi-squared statisticχ2r to test whether the data in several differentfrequency bands are consistent with the matching template [82].Values ofχ2r near unity indicate that the signal is consistent with a coalescence.Ifχ2r is greater than unity,ρðtÞis reweighted asˆρ¼ρ=f½1þðχ2rÞ3 =2g1=6[83,84].The final step enforces coincidence between detectors by selectingevent pairs that occur within a15-ms window and come fromthe same template.The15-ms window is determined by the10-ms intersite propagation time plus5ms for uncertainty inarrival time of weak signals.We rank coincident events basedon the quadrature sumˆρc of theˆρfrom both detectors[45]. To produce background data for this search the SNR maxima of one detector are time shifted and a new set of coincident events is computed.Repeating this procedure ∼107times produces a noise background analysis time equivalent to608000years.To account for the search background noise varying acrossthe target signal space,candidate and background events aredivided into three search classes based on template length.The right panel of Fig.4shows the background for thesearch class of GW150914.The GW150914detection-statistic value ofˆρc¼23.6is larger than any background event,so only an upper bound can be placed on its false alarm rate.Across the three search classes this bound is1in 203000years.This translates to a false alarm probability <2×10−7,corresponding to5.1σ.A second,independent matched-filter analysis that uses adifferent method for estimating the significance of itsevents[85,86],also detected GW150914with identicalsignal parameters and consistent significance.When an event is confidently identified as a real gravitational-wave signal,as for GW150914,the back-ground used to determine the significance of other events is reestimated without the contribution of this event.This is the background distribution shown as a purple line in the right panel of Fig.4.Based on this,the second most significant event has a false alarm rate of1per2.3years and corresponding Poissonian false alarm probability of0.02. Waveform analysis of this event indicates that if it is astrophysical in origin it is also a binary black hole merger[44].VI.SOURCE DISCUSSIONThe matched-filter search is optimized for detecting signals,but it provides only approximate estimates of the source parameters.To refine them we use general relativity-based models[77,78,87,88],some of which include spin precession,and for each model perform a coherent Bayesian analysis to derive posterior distributions of the source parameters[89].The initial and final masses, final spin,distance,and redshift of the source are shown in Table I.The spin of the primary black hole is constrained to be<0.7(90%credible interval)indicating it is not maximally spinning,while the spin of the secondary is only weakly constrained.These source parameters are discussed in detail in[39].The parameter uncertainties include statistical errors and systematic errors from averaging the results of different waveform models.Using the fits to numerical simulations of binary black hole mergers in[92,93],we provide estimates of the mass and spin of the final black hole,the total energy radiated in gravitational waves,and the peak gravitational-wave luminosity[39].The estimated total energy radiated in gravitational waves is3.0þ0.5−0.5M⊙c2.The system reached apeak gravitational-wave luminosity of3.6þ0.5−0.4×1056erg=s,equivalent to200þ30−20M⊙c2=s.Several analyses have been performed to determine whether or not GW150914is consistent with a binary TABLE I.Source parameters for GW150914.We report median values with90%credible intervals that include statistical errors,and systematic errors from averaging the results of different waveform models.Masses are given in the source frame;to convert to the detector frame multiply by(1þz) [90].The source redshift assumes standard cosmology[91]. Primary black hole mass36þ5−4M⊙Secondary black hole mass29þ4−4M⊙Final black hole mass62þ4−4M⊙Final black hole spin0.67þ0.05−0.07 Luminosity distance410þ160−180MpcSource redshift z0.09þ0.03−0.04。

Journal of Machine Learning Research15(2014)3183-3186Submitted6/12;Revised6/13;Published10/14ooDACE Toolbox:A Flexible Object-Oriented Kriging ImplementationIvo Couckuyt∗********************* Tom Dhaene******************* Piet Demeester*********************** Ghent University-iMindsDepartment of Information Technology(INTEC)Gaston Crommenlaan89050Gent,BelgiumEditor:Mikio BraunAbstractWhen analyzing data from computationally expensive simulation codes,surrogate model-ing methods arefirmly established as facilitators for design space exploration,sensitivity analysis,visualization and optimization.Kriging is a popular surrogate modeling tech-nique used for the Design and Analysis of Computer Experiments(DACE).Hence,the past decade Kriging has been the subject of extensive research and many extensions have been proposed,e.g.,co-Kriging,stochastic Kriging,blind Kriging,etc.However,few Krig-ing implementations are publicly available and tailored towards scientists and engineers.Furthermore,no Kriging toolbox exists that unifies several Krigingflavors.This paper addresses this need by presenting an efficient object-oriented Kriging implementation and several Kriging extensions,providing aflexible and easily extendable framework to test and implement new Krigingflavors while reusing as much code as possible.Keywords:Kriging,Gaussian process,co-Kriging,blind Kriging,surrogate modeling, metamodeling,DACE1.IntroductionThis paper is concerned with efficiently solving complex,computational expensive problems using surrogate modeling techniques(Gorissen et al.,2010).Surrogate models,also known as metamodels,are cheap approximation models for computational expensive(black-box) simulations.Surrogate modeling techniques are well-suited to handle,for example,expen-sivefinite element(FE)simulations and computationalfluid dynamic(CFD)simulations.Kriging is a popular surrogate model type to approximate deterministic noise-free data. First conceived by Danie Krige in geostatistics and later introduced for the Design and Analysis of Computer Experiments(DACE)by Sacks et al.(1989),these Gaussian pro-cess(Rasmussen and Williams,2006)based surrogate models are compact and cheap to evaluate,and have proven to be very useful for tasks such as optimization,design space exploration,visualization,prototyping,and sensitivity analysis(Viana et al.,2014).Note ∗.Ivo Couckuyt is a post-doctoral research fellow of FWO-Vlaanderen.Couckuyt,Dhaene and Demeesterthat Kriging surrogate models are primarily known as Gaussian processes in the machine learning community.Except for the utilized terminology there is no difference between the terms and associated methodologies.While Kriging is a popular surrogate model type,not many publicly available,easy-to-use Kriging implementations exist.Many Kriging implementations are outdated and often limited to one specific type of Kriging.Perhaps the most well-known Kriging toolbox is the DACE toolbox1of Lophaven et al.(2002),but,unfortunately,the toolbox has not been updated for some time and only the standard Kriging model is provided.Other freely available Kriging codes include:stochastic Kriging(Staum,2009),2DiceKriging,3 Gaussian processes for Machine Learning(Rasmussen and Nickisch,2010)(GPML),4demo code provided with Forrester et al.(2008),5and the Matlab Krigeage toolbox.6 This paper addresses this need by presenting an object-oriented Kriging implementation and several Kriging extensions,providing aflexible and easily extendable framework to test and implement new Krigingflavors while reusing as much code as possible.2.ooDACE ToolboxThe ooDACE toolbox is an object-oriented Matlab toolbox implementing a variety of Krig-ingflavors and extensions.The most important features and Krigingflavors include:•Simple Kriging,ordinary Kriging,universal Kriging,stochastic Kriging(regression Kriging),blind-and co-Kriging.•Derivatives of the prediction and prediction variance.•Flexible hyperparameter optimization.•Useful utilities include:cross-validation,integrated mean squared error,empirical variogram plot,debug plot of the likelihood surface,robustness-criterion value,etc.•Proper object-oriented design(compatible interface with the DACE toolbox1is avail-able).Documentation of the ooDACE toolbox is provided in the form of a getting started guide (for users),a wiki7and doxygen documentation8(for developers and more advanced users). In addition,the code is well-documented,providing references to research papers where appropriate.A quick-start demo script is provided withfive surrogate modeling use cases, as well as script to run a suite of regression tests.A simplified UML class diagram,showing only the most important public operations, of the toolbox is shown in Figure1.The toolbox is designed with efficiency andflexibil-ity in mind.The process of constructing(and predicting)a Kriging model is decomposed in several smaller,logical steps,e.g.,constructing the correlation matrix,constructing the1.The DACE toolbox can be downloaded at http://www2.imm.dtu.dk/~hbn/dace/.2.The stochastic Kriging toolbox can be downloaded at /.3.The DiceKriging toolbox can be downloaded at /web/packages/DiceKriging/index.html.4.The GPML toolbox can be downloaded at /software/view/263/.5.Demo code of Kriging can be downloaded at //legacy/wileychi/forrester/.6.The Krigeage toolbox can be downloaded at /software/kriging/.7.The wiki documentation of the ooDACE toolbox is found at http://sumowiki.intec.ugent.be/index.php/ooDACE:ooDACE_toolbox.8.The doxygen documentation of the ooDACE toolbox is found at http://sumo.intec.ugent.be/buildbot/ooDACE/doc/.Figure1:Class diagram of the ooDACE toolbox.regression matrix,updating the model,optimizing the parameters,etc.These steps are linked together by higher-level steps,e.g.,fitting the Kriging model and making predic-tions.The basic steps needed for Kriging are implemented as(protected)operations in the BasicGaussianProcess superclass.Implementing a new Kriging type,or extending an existing one,is now done by subclassing the Kriging class of your choice and inheriting the(protected)methods that need to be reimplemented.Similarly,to implement a new hyperparameter optimization strategy it suffices to create a new class inherited from the Optimizer class.To assess the performance of the ooDACE toolbox a comparison between the ooDACE toolbox and the DACE toolbox1is performed using the2D Branin function.To that end,20data sets of increasing size are constructed,each drawn from an uniform random distribution.The number of observations ranges from10to200samples with steps of10 samples.For each data set,a DACE toolbox1model,a ooDACE ordinary Kriging and a ooDACE blind Kriging model have been constructed and the accuracy is measured on a dense test set using the Average Euclidean Error(AEE).Moreover,each test is repeated 1000times to remove any random factor,hence the average accuracy of all repetitions is used.Results are shown in Figure2a.Clearly,the ordinary Kriging model of the ooDACE toolbox consistently outperforms the DACE toolbox for any given sample size,mostly due to a better hyperparameter optimization,while the blind Kriging model is able improve the accuracy even more.3.ApplicationsThe ooDACE Toolbox has already been applied successfully to a wide range of problems, e.g.,optimization of a textile antenna(Couckuyt et al.,2010),identification of the elasticity of the middle-ear drum(Aernouts et al.,2010),etc.In sum,the ooDACE toolbox aims to provide a modern,up to date Kriging framework catered to scientists and age instructions,design documentation,and stable releases can be found at http://sumo.intec.ugent.be/?q=ooDACE.ReferencesJ.Aernouts,I.Couckuyt,K.Crombecq,and J.J.J.Dirckx.Elastic characterization of membranes with a complex shape using point indentation measurements and inverseCouckuyt,Dhaene and Demeester(a)(b)Figure2:(a)Evolution of the average AEE versus the number of samples(Branin function).(b)Landscape plot of the Branin function.modelling.International Journal of Engineering Science,48:599–611,2010.I.Couckuyt,F.Declercq,T.Dhaene,and H.Rogier.Surrogate-based infill optimization applied to electromagnetic problems.Journal of RF and Microwave Computer-Aided Engineering:Advances in Design Optimization of Microwave/RF Circuits and Systems, 20(5):492–501,2010.A.Forrester,A.Sobester,and A.Keane.Engineering Design Via Surrogate Modelling:A Practical Guide.Wiley,Chichester,2008.D.Gorissen,K.Crombecq,I.Couckuyt,P.Demeester,and T.Dhaene.A surrogate modeling and adaptive sampling toolbox for computer based design.Journal of Machine Learning Research,11:2051–2055,2010.URL http://sumo.intec.ugent.be/.S.N.Lophaven,H.B.Nielsen,and J.Søndergaard.Aspects of the Matlab toolbox DACE. Technical report,Informatics and Mathematical Modelling,Technical University of Den-mark,DTU,Richard Petersens Plads,Building321,DK-2800Kgs.Lyngby,2002.C.E.Rasmussen and H.Nickisch.Gaussian processes for machine learning(GPML)toolbox. Journal of Machine Learning Research,11:3011–3015,2010.C.E.Rasmussen and C.K.I.Williams.Gaussian Processes for Machine Learning.MIT Press,2006.J.Sacks,W.J.Welch,T.J.Mitchell,and H.P.Wynn.Design and analysis of computer experiments.Statistical Science,4(4):409–435,1989.J.Staum.Better simulation metamodeling:The why,what,and how of stochastic Kriging. In Proceedings of the Winter Simulation Conference,2009.F.A.C.Viana,T.W.Simpson,V.Balabanov,and V.Toropov.Metamodeling in multi-disciplinary design optimization:How far have we really come?AIAA Journal,52(4): 670–690,2014.。

美国是世界上公司法、证券法研究最为发达的国家之一，在美国法学期刊（Law Review & Journals）上每年发表400多篇以公司法和证券法为主题的论文。

自1994年开始，美国的公司法学者每年会投票从中遴选出10篇左右重要的论文，重印于Corporate Practice Commentator，至2008年，已经评选了15年，计177篇论文入选。

以下是每年入选的论文列表：2008年（以第一作者姓名音序为序）：1.Anabtawi, Iman and Lynn Stout. Fiduciary duties for activist shareholders. 60 Stan. L. Rev. 1255-1308 (2008).2.Brummer, Chris. Corporate law preemption in an age of global capital markets. 81 S. Cal. L. Rev. 1067-1114 (2008).3.Choi, Stephen and Marcel Kahan. The market penalty for mutual fund scandals. 87 B.U. L. Rev. 1021-1057 (2007).4.Choi, Stephen J. and Jill E. Fisch. On beyond CalPERS: Survey evidence on the developing role of public pension funds in corporate governance. 61 V and. L. Rev. 315-354 (2008).5.Cox, James D., Randall S. Thoma s and Lynn Bai. There are plaintiffs and…there are plaintiffs: An empirical analysis of securities class action settlements. 61 V and. L. Rev. 355-386 (2008).6.Henderson, M. Todd. Paying CEOs in bankruptcy: Executive compensation when agency costs are low. 101 Nw. U. L. Rev. 1543-1618 (2007).7.Hu, Henry T.C. and Bernard Black. Equity and debt decoupling and empty voting II: Importance and extensions. 156 U. Pa. L. Rev. 625-739 (2008).8.Kahan, Marcel and Edward Rock. The hanging chads of corporate voting. 96 Geo. L.J. 1227-1281 (2008).9.Strine, Leo E., Jr. Toward common sense and common ground? Reflections on the shared interests of managers and labor in a more rational system of corporate governance. 33 J. Corp. L. 1-20 (2007).10.Subramanian, Guhan. Go-shops vs. no-shops in private equity deals: Evidence and implications.63 Bus. Law. 729-760 (2008).2007年：1.Baker, Tom and Sean J. Griffith. The Missing Monitor in Corporate Governance: The Directors’ & Officers’ Liability Insurer. 95 Geo. L.J. 1795-1842 (2007).2.Bebchuk, Lucian A. The Myth of the Shareholder Franchise. 93 V a. L. Rev. 675-732 (2007).3.Choi, Stephen J. and Robert B. Thompson. Securities Litigation and Its Lawyers: Changes During the First Decade After the PSLRA. 106 Colum. L. Rev. 1489-1533 (2006).4.Coffee, John C., Jr. Reforming the Securities Class Action: An Essay on Deterrence and Its Implementation. 106 Colum. L. Rev. 1534-1586 (2006).5.Cox, James D. and Randall S. Thomas. Does the Plaintiff Matter? An Empirical Analysis of Lead Plaintiffs in Securities Class Actions. 106 Colum. L. Rev. 1587-1640 (2006).6.Eisenberg, Theodore and Geoffrey Miller. Ex Ante Choice of Law and Forum: An Empirical Analysis of Corporate Merger Agreements. 59 V and. L. Rev. 1975-2013 (2006).7.Gordon, Jeffrey N. The Rise of Independent Directors in the United States, 1950-2005: Of Shareholder V alue and Stock Market Prices. 59 Stan. L. Rev. 1465-1568 (2007).8.Kahan, Marcel and Edward B. Rock. Hedge Funds in Corporate Governance and Corporate Control. 155 U. Pa. L. Rev. 1021-1093 (2007).ngevoort, Donald C. The Social Construction of Sarbanes-Oxley. 105 Mich. L. Rev. 1817-1855 (2007).10.Roe, Mark J. Legal Origins, Politics, and Modern Stock Markets. 120 Harv. L. Rev. 460-527 (2006).11.Subramanian, Guhan. Post-Siliconix Freeze-outs: Theory and Evidence. 36 J. Legal Stud. 1-26 (2007). (NOTE: This is an earlier working draft. The published article is not freely available, and at SLW we generally respect the intellectual property rights of others.)2006年：1.Bainbridge, Stephen M. Director Primacy and Shareholder Disempowerment. 119 Harv. L. Rev. 1735-1758 (2006).2.Bebchuk, Lucian A. Letting Shareholders Set the Rules. 119 Harv. L. Rev. 1784-1813 (2006).3.Black, Bernard, Brian Cheffins and Michael Klausner. Outside Director Liability. 58 Stan. L. Rev. 1055-1159 (2006).4.Choi, Stephen J., Jill E. Fisch and A.C. Pritchard. Do Institutions Matter? The Impact of the Lead Plaintiff Provision of the Private Securities Litigation Reform Act. 835.Cox, James D. and Randall S. Thomas. Letting Billions Slip Through Y our Fingers: Empirical Evidence and Legal Implications of the Failure of Financial Institutions to Participate in Securities Class Action Settlements. 58 Stan. L. Rev. 411-454 (2005).6.Gilson, Ronald J. Controlling Shareholders and Corporate Governance: Complicating the Comparative Taxonomy. 119 Harv. L. Rev. 1641-1679 (2006).7.Goshen , Zohar and Gideon Parchomovsky. The Essential Role of Securities Regulation. 55 Duke L.J. 711-782 (2006).8.Hansmann, Henry, Reinier Kraakman and Richard Squire. Law and the Rise of the Firm. 119 Harv. L. Rev. 1333-1403 (2006).9.Hu, Henry T. C. and Bernard Black. Empty V oting and Hidden (Morphable) Ownership: Taxonomy, Implications, and Reforms. 61 Bus. Law. 1011-1070 (2006).10.Kahan, Marcel. The Demand for Corporate Law: Statutory Flexibility, Judicial Quality, or Takeover Protection? 22 J. L. Econ. & Org. 340-365 (2006).11.Kahan, Marcel and Edward Rock. Symbiotic Federalism and the Structure of Corporate Law.58 V and. L. Rev. 1573-1622 (2005).12.Smith, D. Gordon. The Exit Structure of V enture Capital. 53 UCLA L. Rev. 315-356 (2005).2005年：1.Bebchuk, Lucian Arye. The case for increasing shareholder power. 118 Harv. L. Rev. 833-914 (2005).2.Bratton, William W. The new dividend puzzle. 93 Geo. L.J. 845-895 (2005).3.Elhauge, Einer. Sacrificing corporate profits in the public interest. 80 N.Y.U. L. Rev. 733-869 (2005).4.Johnson, . Corporate officers and the business judgment rule. 60 Bus. Law. 439-469 (2005).haupt, Curtis J. In the shadow of Delaware? The rise of hostile takeovers in Japan. 105 Colum. L. Rev. 2171-2216 (2005).6.Ribstein, Larry E. Are partners fiduciaries? 2005 U. Ill. L. Rev. 209-251.7.Roe, Mark J. Delaware?s politics. 118 Harv. L. Rev. 2491-2543 (2005).8.Romano, Roberta. The Sarbanes-Oxley Act and the making of quack corporate governance. 114 Y ale L.J. 1521-1611 (2005).9.Subramanian, Guhan. Fixing freezeouts. 115 Y ale L.J. 2-70 (2005).10.Thompson, Robert B. and Randall S. Thomas. The public and private faces of derivative lawsuits. 57 V and. L. Rev. 1747-1793 (2004).11.Weiss, Elliott J. and J. White. File early, then free ride: How Delaware law (mis)shapes shareholder class actions. 57 V and. L. Rev. 1797-1881 (2004).2004年：1Arlen, Jennifer and Eric Talley. Unregulable defenses and the perils of shareholder choice. 152 U. Pa. L. Rev. 577-666 (2003).2.Bainbridge, Stephen M. The business judgment rule as abstention doctrine. 57 V and. L. Rev. 83-130 (2004).3.Bebchuk, Lucian Arye and Alma Cohen. Firms' decisions where to incorporate. 46 J.L. & Econ. 383-425 (2003).4.Blair, Margaret M. Locking in capital: what corporate law achieved for business organizers in the nineteenth century. 51 UCLA L. Rev. 387-455 (2003).5.Gilson, Ronald J. and Jeffrey N. Gordon. Controlling shareholders. 152 U. Pa. L. Rev. 785-843 (2003).6.Roe, Mark J. Delaware 's competition. 117 Harv. L. Rev. 588-646 (2003).7.Sale, Hillary A. Delaware 's good faith. 89 Cornell L. Rev. 456-495 (2004).8.Stout, Lynn A. The mechanisms of market inefficiency: an introduction to the new finance. 28 J. Corp. L. 635-669 (2003).9.Subramanian, Guhan. Bargaining in the shadow of takeover defenses. 113 Y ale L.J. 621-686 (2003).10.Subramanian, Guhan. The disappearing Delaware effect. 20 J.L. Econ. & Org. 32-59 (2004)11.Thompson, Robert B. and Randall S. Thomas. The new look of shareholder litigation: acquisition-oriented class actions. 57 V and. L. Rev. 133-209 (2004).2003年：1.A yres, Ian and Stephen Choi. Internalizing outsider trading. 101 Mich. L. Rev. 313-408 (2002).2.Bainbridge, Stephen M. Director primacy: The means and ends of corporate governance. 97 Nw. U. L. Rev. 547-606 (2003).3.Bebchuk, Lucian, Alma Cohen and Allen Ferrell. Does the evidence favor state competition in corporate law? 90 Cal. L. Rev. 1775-1821 (2002).4.Bebchuk, Lucian Arye, John C. Coates IV and Guhan Subramanian. The Powerful Antitakeover Force of Staggered Boards: Further findings and a reply to symposium participants. 55 Stan. L. Rev. 885-917 (2002).5.Choi, Stephen J. and Jill E. Fisch. How to fix Wall Street: A voucher financing proposal for securities intermediaries. 113 Y ale L.J. 269-346 (2003).6.Daines, Robert. The incorporation choices of IPO firms. 77 N.Y.U. L. Rev.1559-1611 (2002).7.Gilson, Ronald J. and David M. Schizer. Understanding venture capital structure: A taxexplanation for convertible preferred stock. 116 Harv. L. Rev. 874-916 (2003).8.Kahan, Marcel and Ehud Kamar. The myth of state competition in corporate law. 55 Stan. L. Rev. 679-749 (2002).ngevoort, Donald C. Taming the animal spirits of the stock markets: A behavioral approach to securities regulation. 97 Nw. U. L. Rev. 135-188 (2002).10.Pritchard, A.C. Justice Lewis F. Powell, Jr., and the counterrevolution in the federal securities laws. 52 Duke L.J. 841-949 (2003).11.Thompson, Robert B. and Hillary A. Sale. Securities fraud as corporate governance: Reflections upon federalism. 56 V and. L. Rev. 859-910 (2003).2002年：1.Allen, William T., Jack B. Jacobs and Leo E. Strine, Jr. Function over Form: A Reassessment of Standards of Review in Delaware Corporation Law. 26 Del. J. Corp. L. 859-895 (2001) and 56 Bus. Law. 1287 (2001).2.A yres, Ian and Joe Bankman. Substitutes for Insider Trading. 54 Stan. L. Rev. 235-254 (2001).3.Bebchuk, Lucian Arye, Jesse M. Fried and David I. Walker. Managerial Power and Rent Extraction in the Design of Executive Compensation. 69 U. Chi. L. Rev. 751-846 (2002).4.Bebchuk, Lucian Arye, John C. Coates IV and Guhan Subramanian. The Powerful Antitakeover Force of Staggered Boards: Theory, Evidence, and Policy. 54 Stan. L. Rev. 887-951 (2002).5.Black, Bernard and Reinier Kraakman. Delaware’s Takeover Law: The Uncertain Search for Hidden V alue. 96 Nw. U. L. Rev. 521-566 (2002).6.Bratton, William M. Enron and the Dark Side of Shareholder V alue. 76 Tul. L. Rev. 1275-1361 (2002).7.Coates, John C. IV. Explaining V ariation in Takeover Defenses: Blame the Lawyers. 89 Cal. L. Rev. 1301-1421 (2001).8.Kahan, Marcel and Edward B. Rock. How I Learned to Stop Worrying and Love the Pill: Adaptive Responses to Takeover Law. 69 U. Chi. L. Rev. 871-915 (2002).9.Kahan, Marcel. Rethinking Corporate Bonds: The Trade-off Between Individual and Collective Rights. 77 N.Y.U. L. Rev. 1040-1089 (2002).10.Roe, Mark J. Corporate Law’s Limits. 31 J. Legal Stud. 233-271 (2002).11.Thompson, Robert B. and D. Gordon Smith. Toward a New Theory of the Shareholder Role: "Sacred Space" in Corporate Takeovers. 80 Tex. L. Rev. 261-326 (2001).2001年：1.Black, Bernard S. The legal and institutional preconditions for strong securities markets. 48 UCLA L. Rev. 781-855 (2001).2.Coates, John C. IV. Takeover defenses in the shadow of the pill: a critique of the scientific evidence. 79 Tex. L. Rev. 271-382 (2000).3.Coates, John C. IV and Guhan Subramanian. A buy-side model of M&A lockups: theory and evidence. 53 Stan. L. Rev. 307-396 (2000).4.Coffee, John C., Jr. The rise of dispersed ownership: the roles of law and the state in the separation of ownership and control. 111 Y ale L.J. 1-82 (2001).5.Choi, Stephen J. The unfounded fear of Regulation S: empirical evidence on offshore securities offerings. 50 Duke L.J. 663-751 (2000).6.Daines, Robert and Michael Klausner. Do IPO charters maximize firm value? Antitakeover protection in IPOs. 17 J.L. Econ. & Org. 83-120 (2001).7.Hansmann, Henry and Reinier Kraakman. The essential role of organizational law. 110 Y ale L.J. 387-440 (2000).ngevoort, Donald C. The human nature of corporate boards: law, norms, and the unintended consequences of independence and accountability. 89 Geo. L.J. 797-832 (2001).9.Mahoney, Paul G. The political economy of the Securities Act of 1933. 30 J. Legal Stud. 1-31 (2001).10.Roe, Mark J. Political preconditions to separating ownership from corporate control. 53 Stan. L. Rev. 539-606 (2000).11.Romano, Roberta. Less is more: making institutional investor activism a valuable mechanism of corporate governance. 18 Y ale J. on Reg. 174-251 (2001).2000年：1.Bratton, William W. and Joseph A. McCahery. Comparative Corporate Governance and the Theory of the Firm: The Case Against Global Cross Reference. 38 Colum. J. Transnat’l L. 213-297 (1999).2.Coates, John C. IV. Empirical Evidence on Structural Takeover Defenses: Where Do We Stand?54 U. Miami L. Rev. 783-797 (2000).3.Coffee, John C., Jr. Privatization and Corporate Governance: The Lessons from Securities Market Failure. 25 J. Corp. L. 1-39 (1999).4.Fisch, Jill E. The Peculiar Role of the Delaware Courts in the Competition for Corporate Charters. 68 U. Cin. L. Rev. 1061-1100 (2000).5.Fox, Merritt B. Retained Mandatory Securities Disclosure: Why Issuer Choice Is Not Investor Empowerment. 85 V a. L. Rev. 1335-1419 (1999).6.Fried, Jesse M. Insider Signaling and Insider Trading with Repurchase Tender Offers. 67 U. Chi. L. Rev. 421-477 (2000).7.Gulati, G. Mitu, William A. Klein and Eric M. Zolt. Connected Contracts. 47 UCLA L. Rev. 887-948 (2000).8.Hu, Henry T.C. Faith and Magic: Investor Beliefs and Government Neutrality. 78 Tex. L. Rev. 777-884 (2000).9.Moll, Douglas K. Shareholder Oppression in Close Corporations: The Unanswered Question of Perspective. 53 V and. L. Rev. 749-827 (2000).10.Schizer, David M. Executives and Hedging: The Fragile Legal Foundation of Incentive Compatibility. 100 Colum. L. Rev. 440-504 (2000).11.Smith, Thomas A. The Efficient Norm for Corporate Law: A Neotraditional Interpretation of Fiduciary Duty. 98 Mich. L. Rev. 214-268 (1999).12.Thomas, Randall S. and Kenneth J. Martin. The Determinants of Shareholder V oting on Stock Option Plans. 35 Wake Forest L. Rev. 31-81 (2000).13.Thompson, Robert B. Preemption and Federalism in Corporate Governance: Protecting Shareholder Rights to V ote, Sell, and Sue. 62 Law & Contemp. Probs. 215-242 (1999).1999年（以第一作者姓名音序为序）：1.Bankman, Joseph and Ronald J. Gilson. Why Start-ups? 51 Stan. L. Rev. 289-308 (1999).2.Bhagat, Sanjai and Bernard Black. The Uncertain Relationship Between Board Composition and Firm Performance. 54 Bus. Law. 921-963 (1999).3.Blair, Margaret M. and Lynn A. Stout. A Team Production Theory of Corporate Law. 85 V a. L. Rev. 247-328 (1999).4.Coates, John C., IV. “Fair V alue” As an A voidable Rule of Corporate Law: Minority Discounts in Conflict Transactions. 147 U. Pa. L. Rev. 1251-1359 (1999).5.Coffee, John C., Jr. The Future as History: The Prospects for Global Convergence in Corporate Governance and Its Implications. 93 Nw. U. L. Rev. 641-707 (1999).6.Eisenberg, Melvin A. Corporate Law and Social Norms. 99 Colum. L. Rev. 1253-1292 (1999).7.Hamermesh, Lawrence A. Corporate Democracy and Stockholder-Adopted By-laws: Taking Back the Street? 73 Tul. L. Rev. 409-495 (1998).8.Krawiec, Kimberly D. Derivatives, Corporate Hedging, and Shareholder Wealth: Modigliani-Miller Forty Y ears Later. 1998 U. Ill. L. Rev. 1039-1104.ngevoort, Donald C. Rereading Cady, Roberts: The Ideology and Practice of Insider Trading Regulation. 99 Colum. L. Rev. 1319-1343 (1999).ngevoort, Donald C. Half-Truths: Protecting Mistaken Inferences By Investors and Others.52 Stan. L. Rev. 87-125 (1999).11.Talley, Eric. Turning Servile Opportunities to Gold: A Strategic Analysis of the Corporate Opportunities Doctrine. 108 Y ale L.J. 277-375 (1998).12.Williams, Cynthia A. The Securities and Exchange Commission and Corporate Social Transparency. 112 Harv. L. Rev. 1197-1311 (1999).1998年：1.Carney, William J., The Production of Corporate Law, 71 S. Cal. L. Rev. 715-780 (1998).2.Choi, Stephen, Market Lessons for Gatekeepers, 92 Nw. U. L. Rev. 916-966 (1998).3.Coffee, John C., Jr., Brave New World?: The Impact(s) of the Internet on Modern Securities Regulation. 52 Bus. Law. 1195-1233 (1997).ngevoort, Donald C., Organized Illusions: A Behavioral Theory of Why Corporations Mislead Stock Market Investors (and Cause Other Social Harms). 146 U. Pa. L. Rev. 101-172 (1997).ngevoort, Donald C., The Epistemology of Corporate-Securities Lawyering: Beliefs, Biases and Organizational Behavior. 63 Brook. L. Rev. 629-676 (1997).6.Mann, Ronald J. The Role of Secured Credit in Small-Business Lending. 86 Geo. L.J. 1-44 (1997).haupt, Curtis J., Property Rights in Firms. 84 V a. L. Rev. 1145-1194 (1998).8.Rock, Edward B., Saints and Sinners: How Does Delaware Corporate Law Work? 44 UCLA L. Rev. 1009-1107 (1997).9.Romano, Roberta, Empowering Investors: A Market Approach to Securities Regulation. 107 Y ale L.J. 2359-2430 (1998).10.Schwab, Stewart J. and Randall S. Thomas, Realigning Corporate Governance: Shareholder Activism by Labor Unions. 96 Mich. L. Rev. 1018-1094 (1998).11.Skeel, David A., Jr., An Evolutionary Theory of Corporate Law and Corporate Bankruptcy. 51 V and. L. Rev. 1325-1398 (1998).12.Thomas, Randall S. and Martin, Kenneth J., Should Labor Be Allowed to Make Shareholder Proposals? 73 Wash. L. Rev. 41-80 (1998).1997年：1.Alexander, Janet Cooper, Rethinking Damages in Securities Class Actions, 48 Stan. L. Rev. 1487-1537 (1996).2.Arlen, Jennifer and Kraakman, Reinier, Controlling Corporate Misconduct: An Analysis of Corporate Liability Regimes, 72 N.Y.U. L. Rev. 687-779 (1997).3.Brudney, Victor, Contract and Fiduciary Duty in Corporate Law, 38 B.C. L. Rev. 595-665 (1997).4.Carney, William J., The Political Economy of Competition for Corporate Charters, 26 J. Legal Stud. 303-329 (1997).5.Choi, Stephen J., Company Registration: Toward a Status-Based Antifraud Regime, 64 U. Chi. L. Rev. 567-651 (1997).6.Fox, Merritt B., Securities Disclosure in a Globalizing Market: Who Should Regulate Whom. 95 Mich. L. Rev. 2498-2632 (1997).7.Kahan, Marcel and Klausner, Michael, Lockups and the Market for Corporate Control, 48 Stan. L. Rev. 1539-1571 (1996).8.Mahoney, Paul G., The Exchange as Regulator, 83 V a. L. Rev. 1453-1500 (1997).haupt, Curtis J., The Market for Innovation in the United States and Japan: V enture Capital and the Comparative Corporate Governance Debate, 91 Nw. U.L. Rev. 865-898 (1997).10.Skeel, David A., Jr., The Unanimity Norm in Delaware Corporate Law, 83 V a. L. Rev. 127-175 (1997).1996年：1.Black, Bernard and Reinier Kraakman A Self-Enforcing Model of Corporate Law, 109 Harv. L. Rev. 1911 (1996)2.Gilson, Ronald J. Corporate Governance and Economic Efficiency: When Do Institutions Matter?, 74 Wash. U. L.Q. 327 (1996)3. Hu, Henry T.C. Hedging Expectations: "Derivative Reality" and the Law and Finance of the Corporate Objective, 21 J. Corp. L. 3 (1995)4.Kahan, Marcel & Michael Klausner Path Dependence in Corporate Contracting: Increasing Returns, Herd Behavior and Cognitive Biases, 74 Wash. U. L.Q. 347 (1996)5.Kitch, Edmund W. The Theory and Practice of Securities Disclosure, 61 Brooklyn L. Rev. 763 (1995)ngevoort, Donald C. Selling Hope, Selling Risk: Some Lessons for Law From Behavioral Economics About Stockbrokers and Sophisticated Customers, 84 Cal. L. Rev. 627 (1996)7.Lin, Laura The Effectiveness of Outside Directors as a Corporate Governance Mechanism: Theories and Evidence, 90 Nw. U.L. Rev. 898 (1996)lstein, Ira M. The Professional Board, 50 Bus. Law 1427 (1995)9.Thompson, Robert B. Exit, Liquidity, and Majority Rule: Appraisal's Role in Corporate Law, 84 Geo. L.J. 1 (1995)10.Triantis, George G. and Daniels, Ronald J. The Role of Debt in Interactive Corporate Governance. 83 Cal. L. Rev. 1073 (1995)1995年：公司法：1.Arlen, Jennifer and Deborah M. Weiss A Political Theory of Corporate Taxation,. 105 Y ale L.J. 325-391 (1995).2.Elson, Charles M. The Duty of Care, Compensation, and Stock Ownership, 63 U. Cin. L. Rev. 649 (1995).3.Hu, Henry T.C. Heeding Expectations: "Derivative Reality" and the Law and Finance of the Corporate Objective, 73 Tex. L. Rev. 985-1040 (1995).4.Kahan, Marcel The Qualified Case Against Mandatory Terms in Bonds, 89 Nw. U.L. Rev. 565-622 (1995).5.Klausner, Michael Corporations, Corporate Law, and Networks of Contracts, 81 V a. L. Rev. 757-852 (1995).6.Mitchell, Lawrence E. Cooperation and Constraint in the Modern Corporation: An Inquiry Into the Causes of Corporate Immorality, 73 Tex. L. Rev. 477-537 (1995).7.Siegel, Mary Back to the Future: Appraisal Rights in the Twenty-First Century, 32 Harv. J. on Legis. 79-143 (1995).证券法：1.Grundfest, Joseph A. Why Disimply? 108 Harv. L. Rev. 727-747 (1995).2.Lev, Baruch and Meiring de V illiers Stock Price Crashes and 10b-5 Damages: A Legal Economic, and Policy Analysis, 47 Stan. L. Rev. 7-37 (1994).3.Mahoney, Paul G. Mandatory Disclosure as a Solution to Agency Problems, 62 U. Chi. L. Rev. 1047-1112 (1995).4.Seligman, Joel The Merits Do Matter, 108 Harv. L. Rev. 438 (1994).5.Seligman, Joel The Obsolescence of Wall Street: A Contextual Approach to the Evolving Structure of Federal Securities Regulation, 93 Mich. L. Rev. 649-702 (1995).6.Stout, Lynn A. Are Stock Markets Costly Casinos? Disagreement, Mark Failure, and Securities Regulation, 81 V a. L. Rev. 611 (1995).7.Weiss, Elliott J. and John S. Beckerman Let the Money Do the Monitoring: How Institutional Investors Can Reduce Agency Costs in Securities Class Actions, 104 Y ale L.J. 2053-2127 (1995).1994年：公司法：1.Fraidin, Stephen and Hanson, Jon D. Toward Unlocking Lockups, 103 Y ale L.J. 1739-1834 (1994)2.Gordon, Jeffrey N. Institutions as Relational Investors: A New Look at Cumulative V oting, 94 Colum. L. Rev. 124-192 (1994)3.Karpoff, Jonathan M., and Lott, John R., Jr. The Reputational Penalty Firms Bear From Committing Criminal Fraud, 36 J.L. & Econ. 757-802 (1993)4.Kraakman, Reiner, Park, Hyun, and Shavell, Steven When Are Shareholder Suits in Shareholder Interests?, 82 Geo. L.J. 1733-1775 (1994)5.Mitchell, Lawrence E. Fairness and Trust in Corporate Law, 43 Duke L.J. 425- 491 (1993)6.Oesterle, Dale A. and Palmiter, Alan R. Judicial Schizophrenia in Shareholder V oting Cases, 79 Iowa L. Rev. 485-583 (1994)7. Pound, John The Rise of the Political Model of Corporate Governance and Corporate Control, 68 N.Y.U. L. Rev. 1003-1071 (1993)8.Skeel, David A., Jr. Rethinking the Line Between Corporate Law and Corporate Bankruptcy, 72 Tex. L. Rev. 471-557 (1994)9.Thompson, Robert B. Unpacking Limited Liability: Direct and V icarious Liability of Corporate Participants for Torts of the Enterprise, 47 V and. L. Rev. 1-41 (1994)证券法：1.Alexander, Janet Cooper The V alue of Bad News in Securities Class Actions, 41 UCLA L.Rev. 1421-1469 (1994)2.Bainbridge, Stephen M. Insider Trading Under the Restatement of the Law Governing Lawyers, 19 J. Corp. L. 1-40 (1993)3.Black, Bernard S. and Coffee, John C. Jr. Hail Britannia?: Institutional Investor Behavior Under Limited Regulation, 92 Mich. L. Rev. 1997-2087 (1994)4.Booth, Richard A. The Efficient Market, portfolio Theory, and the Downward Sloping Demand Hypothesis, 68 N.Y.U. L. Rev. 1187-1212 (1993)5.Coffee, John C., Jr. The SEC and the Institutional Investor: A Half-Time Report, 15 Cardozo L. Rev 837-907 (1994)6.Fox, Merritt B. Insider Trading Deterrence V ersus Managerial Incentives: A Unified Theory of Section 16(b), 92 Mich. L. Rev. 2088-2203 (1994)7.Grundfest, Joseph A. Disimplying Private Rights of Action Under the Federal Securities Laws: The Commission's Authority, 107 Harv. L. Rev. 961-1024 (1994)8.Macey, Jonathan R. Administrative Agency Obsolescence and Interest Group Formation: A Case Study of the SEC at Sixty, 15 Cardozo L. Rev. 909-949 (1994)9.Rock, Edward B. Controlling the Dark Side of Relational Investing, 15 Cardozo L. Rev. 987-1031 (1994)。

Design,Optimization and Control of Extractive Distillation for the Separation of Trimethyl Borate −MethanolZhenyu Bao,Weijiang Zhang,Xianbao Cui,and Jiao Xu *Department of Chemical Engineering,Tianjin University,Tianjin 300072,China1.INTRODUCTIONTrimethyl borate (B(OCH 3)3)is an essential reagent for the production of organic boride,sodium borohydride,high-purity boron,antifriction additive of lubricants,stabilizer and plasticizer of polymers,etc.Due to its particular physical and chemical properties,trimethyl borate is also used as a catalyst,dehydrant,auxiliary solvent in brass-welding and high-energy fuel in aircrafts.Natural abundance boron contains 19.8%B-10and 80.2%B-11,1which are the only two stable isotopes.Natural abundance trimethyl borate is normally synthesized by esteri fication of boric acid and methanol.It can also be synthesized by the reaction of borax,sulfuric acid and methanol.The reaction mixtures of the two synthesis methods contain trimethyl borate and methanol,which should be separated for further utilization.As B-10has excellent neutron absorption performance,it is widely used in nuclear power plants,military equipment and medical treatment.High abundance boron used in pressurized water reactors often exists in the form of boric acid.2Production of high B-10abundance boric acid is reported by Han et al.,3which inevitably forms a mixture of trimethyl borate and methanol,and the yield of target product is strongly a ffected by the separation of the mixture.At 101.3kPa,trimethyl borate and methanol form a minimum-boiling azeotrope,the azeotropic point is 328.37K and its azeotropic mixture contains 77−78%(wt %)trimethyl borate.4The azeotrope cannot be separated by ordinary distillation,so special distillation such as extractive distillation and other separation methods are introduced to separate the azeotropic mixture.Sulfuric acid dealcoholization is one of the industrially used methods.The product is an intermediate containing 92%(wt %)trimethyl borate.Then,a salting-out process is employed to get higher purity trimethyl borate (i.e.,98%(wt %)).Usually,lithium chloride is applied.5The separation process consumes a lot of sulfuric acid,which is corrosive,and the salting-out process is inconvenient to carry out.Extractive distillation is an important technique to separate close-boiling point mixture and azeotropes.It is based on the preferential a ffinity of solvent (entrainer)for one or more of the components over the others,which thus alters the relativevolatilities of the feed components.Khoury has given detailed descriptions on such a process.A lot of azeotropes have been separated via extractive distillation;however,up to now,the separation of the azeotrope of trimethyl borate −methanol by extractive distillation has not been reported in open literature.Solvent selection is a major task in extractive distillation.In recent years,many methods have been developed to select proper solvents,such as qualitative judgment,quantitative estimation and experimental determination.Pretel et al.7and Papadopoulos and Linke 8have applied computer-aided molecular design (CAMD)in single solvent design.Moreover,Karunanithi and Achenie 9have studied solvent mixture designs.Dyk and Nieuwoudt 10,11have extended genetic algorithm based CAMD models to design solvents mainly used for the distillation process.Conceptual design and optimization of distillation process have been investigated by Douglas,12Doherty and Malone 13andRodri guez-Donis et al.14Their research e fforts have shown how to select solvents and optimize column sequences.Because dynamic performance proves to be a necessary way to assess the e ffectiveness and controllability of steady state design,Sakizlis 15has presented reviews of integrated design and control methodologies.Shirsat et al.16and Ghaee et al.17have given detailed cases concerning extractive distillation.An e ffective plantwide process control system includes several goals such as safe and smooth process operation and tight control of product quality in front of disturbances.Luyben et al.18−22have presented essential tips and theories on dynamic control using Aspen Plus Dynamics and HYSYS.In their works,development of rigorous simulation of single distillation columns and sequences of columns are concluded;the interaction between steady state design and control is explained;plantwide control with emphasis on selection of control structures for an entire multiunit process is introduced;feedforward,feedback and protective controls areReceived:May 19,2014Revised:August 18,2014Accepted:August 26,2014applied to achieve automatic startup,shutdown and smooth,noninteracting control of column product composition.To fill in the gap of research on trimethyl borate −methanol separation technique,in this study,a proper solvent was selected with a combination of software calculation and experiment determination,and the solvent was comprehensively compared with other potential solvents.Then,the process was optimized through an iteration algorithm.Optimized parameters were further testi fied using the single factor method.At last,control structures were proposed to improve the dynamic stability of the system.2.PROCESS DESIGN AND OPTIMIZATION2.1.Solvent Selection.Solvent selection is the key problem for extractive distillation.The performance of a solvent can be indicated by several indicators based on the activity coe fficient at in finite dilution (γ∞).The commonly used indicators are solvent selectivity,solvent power,relative volatility at in finite dilution in the solvent and performance index.Their de finitions are as follows:9,7βγγ=∞∞solvent selectivity1,S 2,S(1)γ=·∞solvent power SP 1MW MW 2,S 2S(2)αγγ=∞∞P P relative volatility1,21,S 1S2,S2S (3)α=x performance indexPI MW 1,2S ms(4)where x ms is the minimum solvent molar fraction to break theazeotrope.For estimation of γ∞,an activity-coe fficient estimation model called conductor-like screening model for segment activity coe fficient (COSMO-SAC)is used,which is based on the interactions between surface charge distributions of molecules in solution.23Molecules are first transferred from vacuum phase to an ideal conductor,where surface charge is ideally shielded,then transferred to real solvent using polarity factors,rather than directly into condensed phase.COSMO-SAC obtains the molecular segments activity coe fficient using solvation free energy,which avoids the violation of thermodynamic consistency principle in several boundary conditions.24−26Most approaches in selecting solvents are mainly based on the same concept of molecular generation,testing and matching the molecules with speci fied target molecular properties.27−30Considering the distinct polarity di fference between trimethyl borate and methanol,combined with price,toxicity and stabilityconsideration in industrial application,several frequently used strong polar solvents such as dimethyl sulfoxide (DMSO),N ,N -dimethylformamide (DMF),N ,N -dimethylacetamide (DMAC),ethylene glycol (EG),glycerin and N -methyl-2-pyrrolidinone (NMP)were investigated.For each solvent,ten vapor −liquid equilibrium (VLE)data points were experimentally determined using the same method and facilities with Tian et al.32The regression and correlation were made using ChemCAD.The minimum solvent molar fraction (x ms )to break the azeotrope was then obtained,and the result is shown in Table 1.Meanwhile,vapor pressures of trimethyl borate and methanol were calculated using the Antoine equation with coe fficients obtained from literature.33,34As can be seen from Table 1,DMSO has comparatively higher β·SP and PI values than other solvents,which indicates a favorable e ffect on separation of the azeotrope.Moreover,it has a relatively low heat capacity,so it is more energy saving than other solvents.There is an interesting phenomenon that the x ms value is lower for NMP than for DMSO.But with increasing the solvent molar fraction,the solvent e ffect of DMSO increases sharply,which displays a rather high relative volatility when the solvent molar fraction is 0.4(see Figure 1).This obvious tendency was validated by further increase of the DMSO molar fraction.Figure 1shows the XY phase diagram for trimethyl borate −methanol with di fferent solvents (solvent molar fraction is 0.4)at 101.3kPa.Values and curves of EG and glycerin are not shown in Table 1and Figure 1because their miscibility and flowability are bad;moreover,the problem can be aggravated when mixing with trimethyl borate and methanol,as the temperatures of the mixtures are reduced.Table 1.Results of Solvent Selectionsolvent T b (K)γ1,S ∞γ2,S ∞βSP β·SP C P a (J ·mol −1·K −1)x ms α1,2100×PI DMSO 462.15 1.26850.3925 3.2323 1.0449 3.3774149.390.33 1.13 4.37DMF 426.150.98740.5161 1.91330.8494 1.6252146.050.440.87 2.72DMAC 438.150.97060.3478 2.7906 1.0574 2.9508178.20.361.163.71EG 470.35 3.09510.9719 3.18460.5311 1.6913149.8glycerol 563.15 2.89170.9571 3.02130.3635 1.0982218.9NMP476.150.91810.30473.01291.06083.1961412.40.290.95 3.30aData from NIST Chemistry WebBook.31Figure 1.Isobaric VLE diagram for the system of trimethyl borate (1)+methanol (2)+solvent (3)(x 3=0.4)at 101.3kPa.In Figure 1,x 1′represents the mole fraction of trimethyl borate in the liquid phase excluding the solvent,y 1is mole fraction of trimethyl borate in the vapor phase.It can be seen that DMSO can greatly enhance the relative volatility of trimethyl borate and methanol when the molar fraction is 0.4.NMP is another potential solvent with a strong solvent e ffect,especially when separating mixtures with high trimethyl borate concentrations.Whether the mixture of DMSO and NMP can produce better performance needs further investigation.In contrast,DMAC and DMF seem to be not as good as the former two solvents.The overall order of the solvent performance may be resulted from group interaction di fference between the mixture and the N-containing/S-containing compounds.2.2.Feasibility Study.The feasibility study was carried out by Aspen Plus V7.2.The UNIQUAC activity coe fficient model was used to calculate the vapor −liquid equilibrium;however,only the binary-interaction parameters of methanol −DMSO were found in the database.The binary-interaction parameters of trimethyl borate −methanol were obtained from Gmehling ’s handbook,35and the binary-interaction parameters of trimethyl borate −DMSO were correlated from the vapor −liquid equili-brium data calculated by COSMO-SAC.Results are shown in Table 2,and the parameters of UNIQUAC model are 0except for bij and bji.A residue curve map of the ternary system was obtained through Aspen Plus and is shown in Figure 2.As there exists noboundary line in the map,the azeotrope and DMSO is the unstable and stable node,respectively,so the extractive distillation process is feasible.In this figure,black bold lines are the material balance line,F stands for the composition of the material to be separated (fresh feed).F1is the mixture of fresh feed and solvent,which can be separated into D1(overhead of extractive distillation column)and B1(bottom product of extractive distillation column),then B1can be separated into D2(overhead of solvent recovery column)and B2(bottom productof solvent recovery column).D1and D2are trimethyl borate and methanol products,respectively,so it is feasible to separate the azeotrope into pure products with the aid of DMSO.The dashed line across the map is the isovolatility curve.36When DMSO is added,the isovolatility line moves toward the hypotenuse (where the concentration of methanol is zero).Meanwhile,because DMSO lies in the part below the isovolatility curve,where the relative volatility α(C 3H 9BO 3/CH 4O)>1,causing trimethyl borate to go up the column.So,trimethyl borate is the main product of D1.2.3.Partial Optimization.The feed to be separated contains 77%(wt %)trimethyl borate and 23%(wt %)methanol,and the flow rate is 3000kg/h.The product speci fications are the following:D1,99.5%(wt %)trimethyl borate;D2,99.5%(wt %)methanol.The nonoptimized flowsheet is shown in Figure 3.In Figure 3,the variables with interrogation mark behind should be optimized,which will be discussed in the following part.A traditional separation sequence is adopted.As the solvent is recycled for the utmost utilization,a solvent makeup flow is necessarily added to balance the tiny solvent loss.2.3.1.Design of Extractive Column (C1).Because the price of trimethyl borate is much higher than that of methanol,trimethyl borate should be recovered as much as possible.Speci fication for the bottom product is the following:trimethyl borate flow rate should be no more than 0.5kg/h (recovery rate:99.98%).There are five variables in C1that need to be optimized:total stages (N T1);fresh feed stage (N F1);solvent feed stage (N FS );solvent flow rate (S );re flux ratio (RR1).At first,we set the total stages and the feed stage of solvent recovery column at 16and 6,respectively.N T1is given as 62;after trial and error,the minimum energy consumption value (QR1+QR2)is obtained when N F1is set at 48and N FS at 4.Meanwhile,N F1and N FS values other than those always result in disquali fication of products.Similarly,to reach the speci fications,N T1should be no less than 62.That preliminarily con firms the minimum N T1and the proper N F1.Then,concentrations of components in the overhead flow of the extractive column at di fferent S are studied by changing RR1.The results have been plotted in Figure 4,in which the results are obtained when bottom speci fication is satis fied.From Figure 4,we can see that the trimethyl borate concentration reaches a maximum value when RR1=0.4;at the same time,the concentration of impurities is almost the least.Table 2.UNIQUAC Model Parameters of the Ternary Systemcomp,i C 3H 9BO 3C 3H 9BO 3CH 4O comp,j CH 4O DMSO DMSO bij −629.79−162.60129.36bji76.44100.6623.49Figure 2.Residue curve map for trimethyl borate −methanol −DMSO (101.3kPa).Figure 3.Non-optimized flowsheet for the process.Meanwhile,an increase of S results in a higher product purity.However,more solvent also leads to a larger QR1+QR2.The green lines in Figure 4are the purity requirements.As we have de fined that product purity should be 99.5%(wt %),S =24000kg/h is adopted.If we change the purity requirement to 99%(wt %),S can be lowered to 18000kg/h,which can save a lot of energy.Adjusting the operation condition to meet the corresponding purity requirement is suggested when implement-ing the project.The concentration of DMSO is constant when changing S ,because the speci fication is used,but it can be seen that solvent loss reaches the minimum when RR1=0.4.2.3.2.Design of Solvent Recovery Column (C2).Separation of methanol and DMSO takes place in C2.It is relatively easy due to a large relative volatility (76.07,calculated using Aspen Plus built-in data).Methanol is the second target product;99.5%(wt %)of D2is methanol.Most of DMSO enters B2with trace impurity inside (1ppm).There are three variables that need to be optimized:total stages (N T2);feed stage (N F2);re flux ratio (RR2).By fixing the obtained values for C1,we set N T2at 16,then N F2is tested to minimize the energy consumption.Stage 6is selected for N F2,and the initial RR2is estimated by using the same method in the C1design.N T2values other than 16and N F2values other than 6are discussed in the next part.As RR2is manipulated by Aspen Plus to attain the product speci fication,RR2=2.0is used as an initial assignment.2.4.Global Economic Optimization.The tradeo ffbetween economic bene fits,product quality and controllability has been studied by Brengel and Seider,37Luyben and Floudas,38and Palazoglu and Aarkun.39It is a complex job to take every factor into account,so in this study,only utility consumption is considered when calculating operating cost.Total annual cost (TAC)consists of fixed capital investment (FCI)and cost of utilities (C UT ).We will adopt the equation mentioned in Mun o z ’s research:40=+C TAC (10$/a)0.3FCI3UT (eq.5)By observing the proper RR values,we can see that when total stages reduces,FCI drops obviously;however,RR values need to increase to meet the speci fications,which,in turn,raise C UT .TAC is used here to search for the optimum condition.Global optimization is undertaken in the following iterative manner:(1)Fix the pressures at the top of the two columns at 101.3kPa.(2)Give initial estimates for N T2and N F2.(3)Give values for N T1,N F1and N FS .(4)Set S at the minimum value obtained in the partialoptimization part.(5)Change the values of S ,D1and RR1to achieve thespeci fications of C1.(6)Go back to (4)until QR1+QR2is minimized with N T2,N F2,N T1,N F1and N FS fixed.(7)Go back to (3)until TAC is minimized with N T2and N F2fixed.(8)Give values for N T2and N F2.(9)Change the values of D2and RR2to achieve thespeci fications of C2.(10)Go back to (9)until QR1+QR2is minimized with N T1,N F1,N FS ,N T2and N F2fixed.(11)Go back to (4)until TAC is minimized with N T1,N F1andN FS fixed.(12)Go back to (2)until TAC is totally minimized with onlypressures fixed.The optimization procedure is brie fly shown in Figure 5.The results have been listed in Table 3.In Table 3,D 1/2stands for the diameter of the columns.C UT includes the cost of cooling water (0.354$/GJ),steam (14.19$/GJ)and electricity (16.8$/GJ),which are taken from Turton ’s book.41Pumps are considered as electricity-consuming devices.Cost of facilities such as columns,condensers and reboilers are included in FCI,because these large instruments need annual maintenance.In section 2.3(partial optimization),N F1and N FS have been optimized with the minimum energy consumption when N T1=62.Because 62is veri fied to be the optimum value for N T1,then N F1=48and N FS =4are fixed in each case except case 1(N F1=50is optimal).For now,TAC values are calculated to compare the best con figuration of the other variables.Cases 1and 2showFigure 4.Concentration of components in D1at di fferent RR1and S values.how N T1affects the TAC.Herein,N T2and N F2arefixed with the optimum value.As can be seen,TAC is lower when N T1=62, which is also the minimum N T1that can achieve the specifications we set.Likewise,cases3,4and5show how N T2affects TAC.In these cases,the middle part of the columns are selected as N F2.It is clear that N T2=16is superior to N T2=18and14.Cases4,2, and6demonstrate the best N F2should be6.RR1is not the best value we estimate due to surplus purification that we defined.But it does not matter that we choose RR1=0.4,as operation condition,purer products can be obtained in that case.For smaller reflux ratios such as0.33,they can lead to a larger output rate,which is beneficial for a plant.2.5.Analysis of Optimization Result.The single factor method is used to test the validity of the optimization procedure. N T1,N F1,N FS,N T2and N F2are studied besides the above-stated RR and S values.Results are shown in Figure6.In Figure6a,N T1is found to be exactly suitable in achieving the product specification at62,and the purity of trimethyl borate increases with the increase of N T1.Figure6b,c shows that the optimal N F1and N FS values are achieved at48and4,respectively, which are found to be the only values satisfying the given specifications.These locations are influencing because they are either near the bottom or near the top of the column.Increasing N F1or decreasing N FS results in a sharp decrease in product purity.In comparison,a more gradual increase in product purity is observed in C2when N T2is larger than16,as shown in Figure 6d.Moreover,product purity stops increasing when N F2is larger than6,which is demonstrated by Figure6e.The different tendencies between the two columns occur due to the difference in separation difficulties.The more difficult to separate a mixture, the greater impact variables cause on product purity.In conclusion,purity of products can reach the specifications we expected with the optimized values.Diameters of the two columns denote the internal diameter, which is calculated by Aspen inner algorithm and rectified by eq 6:42ρ==F V1max V(6) where V max is the maximum vapor velocity(ft/s)andρV is the corresponding vapor density(lb/ft3).The maximum vapor volumetricflow rate andρV are inquired in hydraulic parameters, then V max is calculated by eq6.The cross-sectional area and internal diameter of the column are thus obtained.By comparing with the column diameter in tray sizing tab,tray spacing is adjusted to make them equal.The tray spacing is0.24m forC1 Figure5.Optimization procedure sequence.Table3.Global Optimization Resultparameters case1case2case3case4case5case6 N T1646262626262N T2161618161416N F1/N FS50/448/448/448/448/448/4 N F2669874RR10.30.330.330.330.330.33RR2 1.73 1.71 1.92 2.32 2.95 1.82D1(m) 1.563 1.565 1.565 1.565 1.565 1.565 D2(m) 1.010 1.011 1.035 1.072 1.137 1.023 Q C1(kW)−261.49−266.30−266.38−266.30−266.30−266.36 Q R1(kW)2449.812448.432448.502448.432448.432448.48 Q C2(kW)−573.41−569.78−614.17−697.23−830.00−593.29 Q R2(kW)805.15807.58859.17935.031060.02831.09C UT(103$)2606.252572.472626.772625.862645.612582.34FCI(103$)747.82723.02737.11729.76724.38724.17 TAC(103$)2830.592789.382847.902844.782862.932799.59and 0.253m for C2.This procedure can ensure the column height is correct,which is essential for economic evaluation.At last,the total solvent loss is calculated by material balance equation,thus the makeup flow rate is 4.077kg/h.A global flowsheet is presented in Figure 7with detailed stream and facility operation information in it.The optimal values of parameters shown in Figure 7are consistent with the case 2values in Table 3.It can be seen that solvent usage is economized,as the S /F ratio is less than 0.14%(wt %)when operating.Even so,the amount of recycled solvent is large,further that utilization of the tremendous energy that the recycled flow possessed is bene ficial.Temperature distribution along the columns and concen-tration of the three components in liquid phase are displayed in Figure 8.Figure 8gives a description on how separation occurs on each stage.For C1,the temperature rises rapidly in the above four stages and turns to be flat between stages 4and 47,then becomes steep again;this is because cold solvent and fresh material are fed on stages 4and 48,respectively.Three sections are dividedapparently and each one has a special function.We may as well label these three sections as s1,s2and s3from column top to bottom.As can be seen from colored lines in Figure 8a,DMSOFigure 6.Correlations between product purity and N T1(a),N F1(b),N FS (c),N T2(d)and N F2(e).Figure 7.Optimal design flowsheet.concentration drops sharply from stage 4to stage 1in s1,as high purity trimethyl borate is required in D1.s2serves to minimize methanol concentration in the overhead,but methanol concentration drops slowly from stage 48to stage 4,which is due to the di fficulty in azeotrope separation.Trimethyl borate is stripped o ffin s3,while methanol and DMSO enter C2.Similarly,but seeming simpler in C2,only a tiny amount of trimethyl borate appears between stages 2and 6,the primary task for this column is to split methanol and DMSO.3.CONTROL SYSTEM DESIGNTo achieve the control purpose,many control strategies with di fferent combination of manipulated variables con figurations have been proposed by Skogestad.43He has elaborated these strategies from theoretical points.More recently,Luyben 42has discussed a variety of control structures with the aid of computer software,one category of them is called “Single end ”control structure where the temperature on one tray is controlled.Moreover,based on Luyben ’s “slope criterion ”42in the determination of control point,temperature measurement is easy,fast and accurate,and composition is closely related with temperature,the stage that has a rapid temperature change (large slope in temperature pro file)can be chosen as an excellent control point.Because this temperature can be used as an obvious re flection of the key component composition variation.Maintaining the temperature of this tray can keep the composition pro file in the column unchanged and thus prevent impurities from entering the withdraw flow.Based on that,stage 61and stage 3are chosen for C1and C2,respectively.Temperature drops on the feed stages are also large,but it is not suitable for control,as the fluctuation of feed conditions can easily destroy the e ffectiveness of the control system.An e ffective and practical control structure needs a precise control point and accurate equipment dimensions.A commonly used heuristic is applied here:42the re flux drum and column base are supposed to provide 10min of liquid holdup,and the length-to-diameter ratio is 2:1.In this case,the re flux drum is 0.72and 0.62m in diameter for C1and C2,respectively.The column base is 1.55and 1.49m in diameter for C1and C2,respectively.3.1.Control Structure with Fixed RRs.As we have obtained RR1and RR2in the optimal operation condition;they are entered in the multipliers.Other manipulated variables and their control actions are listed as follows:(1)Feed flow rate is controlled by adjusting inlet valve(Controller:FC,reverse acting).(2)Solvent flow rate is controlled by adjusting V2B,which is 8times feed flow rate (controller:SFC,direct acting).(3)Temperatures of the certain stage in both columns arecontrolled by adjusting reboiler heat duties (controllers:TC1and TC2,reverse acting).(4)Pressures in both columns are controlled by adjusting thecondenser heat duties (controllers:PC1and PC2,reverse acting).(5)Temperature of recycled solvent is controlled by adjustingheat removal rate of the cooler (controller:TC3,reverse acting).(6)Re flux drum levels in both columns are controlled byadjusting withdraw flow rate of distillates (controllers:LC1D and LC2D,direct acting).(7)Base level in extractive distillation column is controlled byadjusting withdraw flow rate of B1(controller:LC1B,direct acting).(8)Base level in solvent recovery column is controlled byadjusting solvent makeup flow rate (controller:LC2B,reverse acting).In Figure 9,solvent flow rate is manipulated by V2B rather than VM because makeup flow is too small to handle fluctuates when operating.SFC is a cascade controller;it can precisely regulate V2B according to the ratio of SF output (8times feed flow rate)and total solvent flow rate.PID controllers are used here except temperature controllers and pressure controllers.Level controllers are set at K C =2and integral time =9999min,because a proportional-only control is used.Flow controllers use the conventional tuning:K C =0.5,integral time =0.3min.Pressure controllers use default settings.Temperature controllers are PIDincr controllers,and they are tuned with dead time =1min.A Tyreus −Luyben tuning rule is selected to update calculated gains and integral times.Close loop tuning is used with relay amplitude of 5%output range.Detailed tuning parameters are listed in Table 4.Dynamic performance is tested by feed flow rate and feed composition disturbances when time is 0.1min.However,this control structure is not so e ffective to deal with several disturbances.It works well for +20%feed flow rate (3000kg/h →3600kg/h),but it needs 20h to keep the system steady if the feed flow rate decreases 1%(3000kg/h →2970kg/h).For composition disturbances,integrator fails to work when trimethyl borate concentration in the feed changes 1%(both positive and negative).The results for feed flow rate disturbances are shown in Figure 10.It can be seen that variables are brought to new steady state in 3h.New temperatures for both columns are consistent with original ones,but the purity of products is changed.More pronounced is methanol concentration changes in D2whenaFigure 8.Temperature and liquid concentration pro files.。