Integral equations for the spin-weighted spheroidal wave function

a r X i v :0807.2577v 1 [a s t r o -p h ] 16 J u l 2008Rico :An Accurate Cosmological Recombination CodeW. A.Fendt ∗Department of Physics,UIUC,1110W.Green Street,Urbana,IL 61801J.Chluba †Max-Planck-Institut f¨u r Astrophysik,Karl-Schwarzschild-Str.1,85741Garching bei M¨u nchen,GermanyJ.A.Rubi˜n o-Mart´ın ‡Instituto de Astrof´ısica de Canarias (IAC),C/V´ıa L´a ctea,s/n,E-38200,La Laguna,Tenerife (Spain)B. D.Wandelt §Department of Physics,UIUC,1110W.Green Street,Urbana,IL 61801Department of Astronomy,UIUC,1002W.Green Street,Urbana,IL 61801and Center for Advanced Studies,UIUC,912W.Illinois Street,Urbana,IL 61801We present Rico ,a code designed to compute the ionization fraction of the Universe during the epoch of hydrogen and helium recombination with an unprecedented combination of speed and accuracy.This is accomplished by training the machine learning code Pico on the calculations of a multi-level cosmological recombination code which self-consistently includes several physical processes that were neglected previously.After training,Rico is used to fit the free electron fraction as a function of the cosmological parameters.While,for example at low redshifts (z 900),much of the net change in the ionization fraction can be captured by lowering the hydrogen fudge factor in Recfast by about 3%,Rico provides a means of effectively using the accurate ionization history of the full recombination code in the standard cosmological parameter estimation framework without the need to add new or refined fudge factors or functions to a simple recombination model.Within the new approach presented here it is easy to update Rico whenever a more accurate full recombination code becomes available.Once trained,Rico computes the cosmological ionization history with negligible fitting error in ∼10milliseconds,a speed-up of at least 106over the full recombination code that was used here.Also Rico is able to reproduce the ionization history of the full code to a level well below 0.1%,thereby ensuring that the theoretical power spectra of CMB fluctuations can be computed to sufficient accuracy and speed for analysis from upcoming CMB experiments like Planck.Furthermore it will enable cross-checking different recombination codes across cosmological parameter space,a comparison that will be very important in order to assure the accurate interpretation of future cosmic microwave background data.I.INTRODUCTIONPlanck,the third generation of satellite missions to study the Cosmic Microwave Background (CMB),is scheduled for launch later this year.It will make ac-curate measurements of the temperature fluctuations in the CMB out to ℓ∼2500as well as provide a detailed picture of the CMB polarization to ℓ∼1500[1].In order to maximize the knowledge gained about the underlying cosmological parameters from this high resolution experi-ment,it is important that our theoretical calculations are equally precise.Currently,uncertainty in the ionization history of the Universe remains one of the major factors limiting the accuracy of power spectrum calculations [e.g.see 2,3].This fact has recently motived several studies on high precision computations of the cosmological hydrogenand2levels of the hydrogen and helium atoms.As an accurate calculation of the ionization history of a single cosmo-logical model using the current version of our full multi-level recombination code takes many hours or even days of computational time,standard parameter estimation methods require some approximation of the ionization history that is fast to evaluate.Currently the standard method used to evaluate the ionization fraction is the Recfast code[4],which models the hydrogen and he-lium atoms as effective3-level systems.In order to ap-proximate the computations done using their multi-level recombination code[5],the authors introduced a fudge factor,f H,to artificially speed up hydrogen recombina-tion by about14%.More recently,Wong et al.[6]up-dated Recfast to include a second fudge factor,b He,to modify the recombination of helium.This change,which allows speeding up He i recombination,was motivated by the results of Switzer and Hirata[7],Kholupenko et al.[8]and Rubi˜n o-Mart´ın et al.[9].It is intended to include the effect of the hydrogen continuum opacity during neu-tral helium recombination(see Sect.II B2).Here we provide a different approach to computing a fast and accurate approximation of the recombination history.Our method uses a regression code based on Pico[10,11],which we will henceforth refer to as Rico, to model the ionization history,X e(z)=N e(z)/N H(z), as a function of the cosmological parameters.Here N e denotes the number density of free electrons,and N H is the number density of hydrogen nuclei in our Uni-verse.Since X e is a smooth function,the polynomial approximation used by Rico is extremely accurate,and one retains nearly all of the information contained in the multi-level code used for training,the training code.As the accuracy of the training code is the major de-terminant of the accuracy of Rico’s approximation,in Sect.II we review the physical processes included in our training code as well as an assessment of the processes that must still be considered and are currently under in-vestigation.In Sect.III we compare the ionization his-tory and corresponding power spectra computed with the multi-level recombination code used for the training of Rico to the different versions of Recfast.Also we in-clude a short discussion of how one could modify the ap-proach employed by Recfast to capture the new physics included in our full code.In Sect.IV we show that the simple regression scheme used by Rico can accurately model the current full multi-level recombination calcula-tion,giving power spectra that are sufficiently accurate for even a cosmic variance limited experiment. Although the training code discussed here does take into account most of the important corrections discussed in the literature so far,it does not yet solve the prob-lem completely.Still more physical processes should be included and the results validated by independent codes. As Rico can be trained equally well on any recombina-tion code it can facilitate crosschecks between those codes as well as a study of how approximations in these codes propagate to cosmological parameter constraints.z∆N e/N ein%FIG.1:Corrections to the hydrogen recombination history.For the curves labeled’l-splitting’and’our most recent 100-shell computation’we compared with the Recfast-code v1.3after removing all of the switches(see Sect.III B).In the other cases we compared with the solution obtained with our code only including the non-equilibrium effects.Our target is to ensure that the ionization history can ultimately be calculated to sufficient precision to avoid biasing parameter estimation from the next generation of CMB experiments.By making it easy to propagate ad-vances in the calculation of the ionization history through to predictions of the CMB power spectra with Rico,fu-ture development can focus on the physics of recombina-tion and to a lesser degree on the computational efficiency of the training code.Rico and its future updates will be made available at /rico.II.DESCRIPTION OF THE MULTI-LEVEL RECOMBINATION CODE USED FORTRAININGIn this section we discuss the previously neglected physical processes that are accounted for in the multi-level recombination code1which here is used to train Rico.As we will show,the addition of these processes lead to small but significant changes in the CMB power spectra(see Sect.III).A comparison of the impact of all the included processes on the recombination history is illustrated in Fig.1and2.In Sect.II C we also give3 an overview and short discussion of processes that willbe included in a future update of Rico.A.Processes included during Hydrogenrecombination1.Non-equilibrium populations of angular momentumsub-statesWhile the previous work of Seager et al.[5]evolved theindividual energy levels of Hydrogen,thereby allowingdepartures from Saha-equilibrium,it was assumed thatfor principal quantum number n>2the angular momen-tum sub-states were in full statistical equilibrium witheach other.Generally,equilibrium between l sub-statesis maintained through collisional interactions.Duringrecombination,however,collisions are much weaker thanradiative processes,leading to a departure from full sta-tistical equilibrium[12,13].In the computations of Chluba et al.[13]this lead toa∼5−10%increase in the free electron fraction,X e,at very low redshift(z 400).Furthermore,a∼0.6%decrease in X e at decoupling was found in comparisonwith Recfast.As pointed out in their work this lattercorrection was not connected with the departures fromfull statistical equilibrium but has been interpreted as ahint towards limitations of the effective3-level approachused in Recfast[13].For hydrogen here we follow the evolution of up to100l-resolved shells.We also include the l-changing col-lisional rates,however their impact on the recombina-tion history is very small( 0.01%with a maximum atz∼900).This problem involves the simultaneous evo-lution of up to∼5000very stiffdifferential equations;atask that is numerically very cumbersome even on mod-ern machines.The detailed treatment of this process isone of the most difficult and time-consuming aspects ofthe whole training code.Here we have now also updated all the physical con-stants according to the Nist-database2.We then foundthat the0.6%decrease in X e at decoupling practicallydisappeared.This is mainly because previously the ef-fect of the reduced mass of hydrogen,which leads to a∼0.05%correction of the ionization potential,was ne-glected in the code of Chluba et al.[13].We have also improved the evaluation of the photoion-ization and recombination rates,which for many levelstakes rather long.In Chluba et al.[13]the recombina-tion rates for all levels were tabulated before the actualcomputation and detailed balance was used to infer thephotoionization rates each time the system is evaluated.This treatment is possible as long as the photon andelectron temperatures do not depart significantly from3About57%of all hydrogen atoms became neutral via this route,and only∼43%through the Lyman-αline[17]4CMB photons in the Universe leads to stimulated transi-tions of the2s-state to the ground level.Assuming that the CMB is given by a pure blackbody with temperature Tγ=2.725(1+z)K,it was shown that the induced emis-sionleads to a∼1%increase in the2s→1s transition rate[18].This increase has the effect of speeding up hy-drogen recombination leading to a maximum change in the free electron fraction of about1.3%at z∼1050(see Fig.1,curve labeled’induced2s-decay’).Since only the pure CMB blackbody contributes signif-icantly to the stimulated2s→1s transition,it is possible to tabulate the effective transition rate as a function of temperature before the computation.However,for the inverse process(see Sect.II A3)the Lyman-αspectral distortion is very important,so we also include it here for consistency.This is accomplished byfirst running a5-shell computation of the hydrogen recombination to obtain a sufficient estimate of the Lyman-αdistortionwithin the considered cosmology.The effective2s→1s transition rate is then tabulated at5000redshift points and log-interpolated.The integrals were evaluated with relative accuracyǫ∼10−8.3.Feedback of Lyman-αphotons on the effective1s→2sabsorption rate of HydrogenDue to the large number of photons produced in Lyman-αtransitions during hydrogen recombination there is a huge excess of radiation over the background spectrum in the far Wien tail of the CMB.As noted by Kholupenko and Ivanchik[19],after some redshifting, these super-Planckian photons are able to re-excite elec-trons from the ground-state to the2s-level,and therefore increase the effective1s→2s absorption rate.This de-crease in the rate of hydrogen recombination(see Fig.1, curve labeled’1s-2s Lyman-αfeedback’)practically can-cels the effect of stimulated2s-decays(see Sect.II A2), leading to a net∼0.6%increase in the free electron frac-tion at z∼980(see Fig.1,curve labeled’induced2s-decay+1s-2s Lyman-αfeedback’).While the low red-shift behavior of this last result differs slightly from the curve given in Kholupenko and Ivanchik[19],it was also recently obtained by Hirata[20].As explained in Sect.II A2,in order to include the Lyman-αfeedback we run a5-shell computation of hy-drogen recombination to obtain a sufficient estimate of the Lyman-αdistortion and then tabulate the1s→2s absorption rate.We checked that using more shells for the simple run does not affect the results significantly.4.Feedback within the H i Lyman-seriesDue to redshifting,all of the Lyman-series photons emitted in the transition of electrons from levels with n>2have to pass through the next lower-lying Ly-man transition,leading to additional feedback correc-zN e/NHFIG.2:Corrections to helium recombination history.For the inlay plot we compared the ionization history with the standard Recfast-code v1.3after removing all of theswitches(see Sect.III B).tions like in the case of Lyman-αabsorption in the2s-1s two-photon continuum[21].However,in this case the photons connected with Ly n are completely absorbed by the Ly(n−1)resonance and eventually all Ly n photons are converted into Lyman-αor2s-1s two-photon decay quanta.This process delays hydrogen recombination, leading to a maximal correction to the ionization history of∆N e/N e∼0.22%at z∼1050(see Fig.1,curve la-beled’Lyman-series feedback’).This result was recently confirmed by Kholupenko et al.[22].To include this feedback we save the spectral distor-tion due to emission of photons in the Lyman-series up to some given n feed(typically we use n feed∼5)after each time step.We then include the distortion in the evaluation of the excitation rate of the next lower-lying Lyman-series transition,assuming that all photons en-tering the resonance are absorbed.Due to the huge op-tical depth in the Lyman-series this procedure is well-justified.In this way the Lyman-series feedback always works like Ly n→Ly(n−1).Note that due to the dis-tance between neighboring Lyman-series resonances,the feedback occurs after the time it takes to redshift into the next transition.B.Processes included during HeliumrecombinationFor He ii-recombination a list of processes similar to hydrogen can be formulated.In particular C.Hirata and E.Switzer have rigorously studied these and several ad-ditional processes[7,23,24].However,since the details5of helium recombination history are not strongly propa-gating to the computation of the CMB power spectra,as afirst step here we only include the two dominant cor-rections,which actually are due to processes that have no analog in the case of hydrogen recombination.Ad-ditional corrections(see discussion in Sect.II C)will be taken into account for a later release of Rico.One should also mention that for the computation in-volving helium we only include5j-resolved shells for the singlet and triplet states[for more detail see9].We found that within the current numerical approach and related approximations this leads to sufficient precision (see Sect.II C1).1.Spin-forbidden He i-23P1−11S0transitionLike for the Lyman-αphotons in hydrogen,due to the low expansion rate of the Universe,He i-21P1−11S0pho-tons have a very hard time escaping from the interac-tion with neutral helium during He ii→He i recombi-nation.Therefore also the recombination of helium is strongly delayed as compared to the Saha case,and the He i-21S0−11S0two-photon transition[A2γ∼51.4s−1; 25]plays a very important role in defining the rate of He i-recombination.Within the standard computation the He i-21S0−11S0two-photon decay channel allowed ∼31%of helium atoms to recombine,while the He i-21P1−11S0transition contributed∼69%[e.g.see table in26].Dubrovich and Grachev[27]realized that in addition to the He i-21S0−11S0two-photon transition,the spin-forbidden He i-23P1−11S0transition helps helium to re-combine even though it also has a very low probability [A∼177s−1;28,29]as compared to the He i-21P1−11S0 resonance(A∼1.80×109s−1).This process,which was not included in the original version of Recfast[4], leads to a∼1.1%decrease in the free electron fraction at z∼1750.Note that in terms of the helium recombi-nation history this corresponds to a very large effect(see Fig.2).By including the spin-forbidden He i-23P1−11S0tran-sitions into the recombination code onefinds that about ∼17%of all helium atom became neutral through the He i-21S0−11S0two-photon transition,∼40%via the He i-21P1−11S0transition and∼43%through this new channel[26].The comparison shows that the spin-forbidden He i-23P1−11S0transition is actually one of the dominant channels for helium recombination.This conclusion is true even when the hydrogen continuum opacity-effect(see section II B2)is also included.In this case only8%of all helium atoms recombine via the He i-21S0−11S0two-photon decay,while the rest is shared between the He i-21P1−11S0and spin-forbidden He i-23P1−11S0transition,with a∼2%contribution from other direct transitions to the ground-state[9].This process can be self-consistently included by sim-ply adding another term similar to the normal resonance transitions to the set of rate equations[27],a modifica-tion that is now also accounted for in Recfast v1.4[6]. Note that for this transition the Sobolev-escape proba-bility should be taken into account,since the Sobolev optical depth still reaches∼3close to its maximum.2.Absorption of He i photons by neutral hydrogenIn addition,the recombination of neutral helium is sped up due to the absorption of He i-21P1−11S0and 23P1−11S0-photons by the tiny fraction of neutral hy-drogen that is already present at redshifts z 2400. While photons are emitted close to optically thick He i-resonances,a small part of them can be absorbed in the Lyman-continuum of hydrogen,allowing additional elec-trons to settle to the ground state of helium.This pro-cess was mentioned by P.J.E.Peebles in the mid90’s [see remark in2],but only recently taken into account by Switzer and Hirata[23]and others[8,9].It leads to the largest correction during the epoch of He ii→He i-recombination that has been investigated so far in the literature,strongly accelerating the recombination of he-lium below z∼2000(see Fig.2,thick solid line). Since this problem is connected with important de-tails in the radiative transfer of optically thick He i-21P1−11S0and23P1−11S0,it is a very time-consuming task to solve self-consistently.Switzer and Hirata[23] studied this problem using a Monte-Carlo approach to solve the quasi-stationary line-transfer problem.They tabulated the escape probability for different values of the Sobolev-optical depth and hydrogen continuum opac-ity and then interpolated these in the actual recombina-tion calculation.Independently,Kholupenko et al.[8] investigated this problem using a simplified analytical approach,which slightly underestimated the acceleration of helium recombination.They provided simple expres-sions for the correction to the Sobolev escape probability of the He i-21P1−11S0and23P1−11S0,which were later also used by Wong et al.[6]to fudge the helium re-combination history.Rubi˜n o-Mart´ın et al.[9]confirmed the results of Switzer and Hirata[23]using the results of a Fokker-Planck approach to solve the line-transfer problem[30].They also gave a simple1D-integral which neglects the redistribution of photons by resonance scat-tering but still reproduces the correction to the escape probability rather well.Very recently Kholupenko et al.[31]reconsidered this problem in more detail analytically and found good agreement with Switzer and Hirata[23] and Rubi˜n o-Mart´ın et al.[9].In this work we follow the approach of Rubi˜n o-Mart´ın et al.[9]and include the acceleration in the helium recom-bination history using their1D-integral expression[see Eq.(B.3)in9]for the correction to the escape probabil-ity.Here we do not include the redshift dependent fudge-function that was introduced by Rubi˜n o-Mart´ın et al.[9] to account for the additional acceleration because of par-tial redistribution.In the future we plan to solve the full6time-dependent problem including the radiative transferin the optically thick lines,a task that will also become necessary for hydrogen recombination.C.Additional processes that will be considered inthe futureThere are a number of additional processes,both forhydrogen and helium recombination,that have already been addressed in the literature,but were not takeninto account in the training code.As mentioned above, C.Hirata and E.Switzer have rigorously studied severaladditional subtle processes(e.g.feedback;two-photon processes;effect of electron scattering;isotope shift of 3He)within the context of helium recombination[see 7,23,24].Also for hydrogen recombination additionalprocesses have been discussed,however not all of these studies are concluded yet.Some of these additional pro-cesses lead to rather small additional corrections,but others may still be important.Below we give a short overview about some of the work that must still be done. We are planning to make a careful survey in the near future.1.Incompleteness of the atomic modelsAt low redshifts(z 800)the rate of hydrogen recom-bination is strongly controlled by the effective recombina-tion coefficient,which itself depends on the completeness of the atomic model of hydrogen.Here we included only up to100shells in our computations,but as pointed out earlier[13],in terms of the ionization history this may still not be enough,and extension to∼200−300shells may be required.However,for the computation of the CMB power spectra these corrections probably are not very important,as even with75shells rather converged results seem to be obtained(see Sect.IV).Pushing to a larger number of l-resolved shells is not trivial,and also more accurate collisional rates may be required.As mentioned above,for the computation involving he-lium we only include5j-resolved shells.This approxima-tion seems to be sufficient at the current level of precision since,unlike the case of hydrogen recombination,the dy-namics of helium recombination are much less controlled by the effective recombination rate,which is strongly con-nected with the completeness of the atomic model and can be fudged to some level.For helium recombination the escape of photons is much more crucial.2.Escape of H i Lyman-αphotonsThe escape of photons from the optically thick H i Lyman-αline is usually modeled using the Sobolev-approximation.The validity of this approximation dur-ing hydrogen recombination has been investigated several times[e.g.32,33],but at the one percent-level a full con-firmation is still necessary.For example,Grachev and Dubrovich[34]recently claimed that the line-recoil effect leads to a∼1.3%speed-up of hydrogen recombination. Within the quasi-stationary approach this results seems robust[30].Given the importance of this problem to the calculation of the power-spectra,independent checks are necessary.This is now being investigated in detail[30].3.Two-photon transitions from higher levels Along with the2s→1s two-photon transition there are also allowed two-photon transitions from higher levels to the ground state.Within the context of hydrogen and helium recombination these transitions werefirst studied by Dubrovich and Grachev[27].They predicted a∼5%decrease in the free electron fraction at z∼1200, however in the computations of the effective two-photon decay rates of the n s and n d-levels they only included thefirst non-resonant term(i.e.due to the dipole matrix element connecting n s/n d→n p)in the infinite sum over intermediate states.Using rate coefficients for the two-photon decay of the 3s and3d-levels in hydrogen as computed by Cresser et al.[35],Wong and Scott[26]found that Dubrovich and Grachev[27]overestimated the impact of two-photon transitions on the ionization history by about one order of magnitude.However,the calculation of Cresser et al.[35]was incomplete.For example they did not include the largest non-resonant term(due to the dipole matrix element connecting n s/n d→n p)in their calculations [36].Also physically it is very difficult,if not impossible, to separate the‘pure’two-photon decay rate from the resonant contributions[see discussions in20,36,37],e.g. because of non-classical interference effects.Later this problem was reinvestigated in more detail, and a lower limit for the impact of two-photon decays during hydrogen recombination was derived,implying that a decrease of slightly more than0.3−0.5%in the free electron fraction at z∼1150can still be expected [36].This estimate was obtained by taking into account departures of the full two-photon line profiles from the Lorentzian shape in the very distant,optically thin part of the red wing of the Lyman-αline.According to these computations,the two-photon decays from s-states seem to slow hydrogen recombination down,while those from d-states speed it up.In addition it was shown that the slight net acceleration of hydrogen recombination seems to be dominated by the3s and3d contribution[36]. However,it was pointed out that thefinal answer can be only given using a full radiative transfer computation, which also takes into account the effects of partial fre-quency redistribution and the feedback of photons from the blue side of the Lyman-αresonance[e.g.see36]. Very recently Hirata[20]showed that including these as-pects of the problem,along with the induced2s-decay and the feedback of Lyman-αphotons on the effective7 1s→2s absorption rate,modifications in the ionizationhistory of the order of±1.3%can be expected.Also2s→1s Raman-scattering seems to play some role.Thiswould be a very important effect and we are currentlyinvestigating it.Also during helium recombination two-photon transi-tions and Raman-scattering are important.Hirata andSwitzer[24]investigated the effect of these,as well as thestimulated21S0→11S0two-photon decay,feedback byphotons from the He i-21P1−11S0and spin-forbiddenHe i-23P1−11S0transition on the11S0→21S0absorp-tion rate.However,the overall changes due to the two-photon transitions and Raman-scattering seem to be in-significant[24].4.Feedback due to helium linesThe feedback of high frequency photons released dur-ing helium recombination should also affect the dynam-ics of hydrogen recombination.Here it is interesting that most of the high energy photons from He iii→He iiwill be reprocessed by neutral helium before they can directly affect hydrogen.Since the number of high fre-quencyγ(He ii)photons is comparable with the number of helium atoms this should still have a rather strong ef-fect.Note that although the ionization history during He iii→He ii-recombination has a very small impact on the CMB power-spectra,the exact time-dependenceof theγ(He ii)→He i feedback is related to how fast He iii→He ii-recombination occurred,so that it may still deserve careful investigation.Also feedback of the He i line during helium recombi-nation still has to be included,and was shown to have a notable delaying effect[23].However in particular thefeedback of He i-21P1−11S0photons on the23P1−11S0 transition may in addition require a fully time-dependenttreatment,since both resonances are only separated by ∼1%in frequency,or∼600Doppler widths relative to the He i-21P1−11S0resonance at z∼2500.Finally,those photons emitted by neutral helium candirectly feedback on hydrogen,but in order to take this feedback into account more detailed computations of the helium recombination spectrum are required.This pro-cess could still affect the hydrogen recombination history on a level exceeding0.1%.5.Details of the primordial chemistry at low redshiftsFor the initial computations in connection with Rec-fast[5]some aspects of the primordial chemistry were included.These have not yet been taken into account in our code,but are expected to have only a small effect, mainly at redshifts z 200.However,it is in principle easy to include these,as well as extending the chemical network using updated rate coefficients[e.g.see38].-3-2-11230 500 1000 1500 2000 2500 3000∆Cl/C lin%multipole lRecfast v1.3Recfast v1.4Temperature (TT)Polarization (EE)FIG.3:Comparison of the angular power spectra computed using the ionization history of our most recent100shell code versus the default settings in Recfast v1.3and v1.4.The two black lines denote our performance benchmark of±3/ℓwhich is roughly2times cosmic variance for the special case offitting a single amplitude parameter(see also Seljak et al.[3]).6.Other small correctionsSimilar to the effect of3He on helium recombination [7],also deuterium should affect hydrogen recombination. Since due to the isotope shift the deuterium Lyman-αres-onance is on the blue side of the hydrogen Lyman-αline, this should slow down hydrogen recombination.How-ever,since the abundance of deuterium is so small,and because at z∼1100the shift is only∼12Doppler-width of the Lyman-αline,this process is probably negligible. For the feedback of the CMB spectral distortion gen-erated during recombination on the1s→2s-absorption rate we only took into account the hydrogen Lyman-αdistortion.However,the high frequency part of the2s-distortion itself may also have to be included,probably leading to another small correction to the ionization his-tory.PARISON WITH RECF ASTThe net effect of inclusion of the physics listed in Sects.II A and II B leads to an increase in the free elec-tron fraction during hydrogen recombination.Close to the maximum of the Thompson visibility function,which is most relevant for computations of the CMB power spectra,one can observe a0.7%correction.At low z the error increases to∼3%at z∼200.During he-lium recombination the correction reaches about−3%at z∼1800(see Fig.2for details).In Figure3we show the corresponding change in the temperature and polarization power spectra when com-paring with the output of Recfast v1.3and v1.4us-ing the default settings.We used the publicly avail-。

Journal of Engineering Mechanics ASCE,Vol.134(2008),in press What Did and Did not Cause Collapseof WTC Twin Towers in New YorkZdenˇe k P.Baˇz ant,Jia-Liang Le,Frank R.Greening and David B.Benson Structural Engineering Report No.07-05/C605cDepartment of Civil and Environmental EngineeringNorthwestern UniversityEvanston,Illinois60208,USAMay27,2007Revised June22,December15,2007,and March31,2008What Did and Did Not Cause Collapseof WTC Twin Towers in New York Zdenˇe k P.Baˇz ant1,Hon.M.ASCE,Jia-Liang Le2,Frank R.Greening3,and David B.Benson4 Abstract:Previous analysis of progressive collapse showed that gravity alone suffices to explain the overall collapse of the World Trade Center(WTC)towers.However,it remains to be checked whether the recent allegations of controlled demolition have any scientific merit.The present analysis proves that they do not.The video record available for thefirst few seconds of collapse is shown to agree with the motion history calculated from the differential equation of progressive collapse but, despite uncertain values of some parameters,it is totally out of range of the free fall hypothesis,on which these allegations rest.It is shown that the observed size range(0.01mm—0.1mm)of the dust particles of pulverized concrete is consistent with the theory of comminution caused by impact, and that less than10%of the total gravitational energy,converted to kinetic energy,sufficed to produce this dust(whereas more than150tons of TNT per tower would have to be installed,into many small holes drilled into concrete,to produce the same pulverization).The air ejected from the building by gravitational collapse must have attained,near the ground,the speed of almost500mph (or223m/s,or803km/h)on the average,andfluctuations must have reached the speed of sound. This explains the loud booms and wide spreading of pulverized concrete and other fragments,and shows that the lower margin of the dust cloud could not have coincided with the crushing front.The resisting upward forces due to pulverization and to ejection of air,dust and solid fragments,neglected in previous studies,are found to be indeed negligible during thefirst few seconds of collapse but not insignificant near the end of crush-down.The calculated crush-down duration is found to match a logical interpretation of seismic record,while the free fall duration grossly disagrees with this record.IntroductionTo structural engineers,the collapse of the World Trade Center(WTC)towers on9/11/2001 came as the greatest surprise since the collapse of Tacoma Narrows Bridge in1940.Immediately after the aircraft impact,the structural frame behaved as expected,but not after thefire.To explain the collapse,it was proposed(on September13,2001;Baˇz ant2001;Baˇz ant and Zhou2002)that viscoplastic buckling of heated and overloaded columns caused the top part of tower to fall through the height of at least one story,and then shown that the kinetic energy of the impact on the lower part must have exceeded the energy absorption capacity of the lower part by an order of magnitude.A meticulous investigation of unprecedented scope and detail, conducted by S.Shyam Sunder’s team at the National Institute of Standards and Technology (NIST2005),supports this explanation.Although NIST did not analyze the overall process of dynamic progressive collapse below thefire zone,it verified a sequence of effects that triggered the collapse:(1)scraping of much of steel insulation byflying objects during aircraft impact (without which the towers would not have collapsed,as concluded by NIST);(2)cutting of many columns,and damage with large deflections of others during aircraft impact;(3)subsequent load redistributions among columns;(4)sagging of heatedfloor trusses and their catenary action, evidenced by multistory inward bowing of perimeter columns;and(5)viscoplastic buckling of heated,damaged and overloaded columns.Universally though has the foregoing explanation of collapse been accepted by the commu-nities of structural engineers and structural mechanics researchers,some outside critics have nevertheless exploited various unexplained observations to disseminate allegations of controlled demolition.The objective of this paper,based on the report by Baˇz ant et al.(2007),is to ex-amine whether those allegations might be scientifically justifiable,and to show that the concept of gravity-driven collapse does not conflict with any observations.1McCormick Institute Professor and W.P.Murphy Professor of Civil Engineering and Materi-als Science,Northwestern University,2145Sheridan Road,CEE/A135,Evanston,Illinois60208;z-bazant@.2Graduate Research Assistant,Northwestern University.3Engineering Consultant,Hamilton,Ontario L8S3X7.4Professor Emeritus,School of Electrical Engrg.and Computer Science,Washington State University,Pull-man,WA99164.NIST examined paint cracking on16perimeter columns and only3of them showed evidence of temperature>250◦C.From annealing studies,the microstructures of recovered steel samples show no evidence of exposure to temperature above600◦C for any significant time(NIST2005, part NCSTAR1-3,p.132)(i.e.,>15min.).Nevertheless,evidence of very high temperature (around800◦C)was found in one unidentified steel column,though it is unlikely that this column was located in thefire stories and experienced high temperature prior to the collapse (NIST2005,part NCSTAR1-3C p.229).Note,though,that only1%of the columns from the fire stories were examined.Consequently,NIST cautioned that thefindings from paint cracking test and annealing studies are not indicative of the steel temperature in thefire stories.Thus, although very high steel temperature are likely,there is no direct evidence.But are high steel temperatures really necessary to explain collapse?Not really.The initial speculation that very high temperatures were necessary to explain collapse must be now revised since tests revealed a strong temperature effect on the yield strength of the steel used.The tests by NIST(2005,part NCSTAR1-3D,p.135,Fig.6-6) showed that,at temperatures150◦C,250◦C and350◦C,the yield strength of the steel used in thefire stories decreased by12%,19%and25%,respectively.These reductions apply to normal durations of laboratory strength tests(up to several minutes).Since the thermally activated decrease of yield stress is a time-dependent process,the yield strength decrease must have been even greater for the heating durations in the towers,which were of the order of one hour.These effects of heating are further documented by the recentfire tests of Zeng et al. (2003),which showed that structural steel columns under a sustained load of50%to70%of their cold strength collapse when heated to250◦C.Although a detailed computer analysis of columns stresses after aircraft impact is certainly possible,it would be quite tedious and demanding,and has not been carried out by NIST. Nevertheless,it can easily be explained that the stress in some surviving columns most likely exceeded88%of their cold strengthσ0.In that case,any steel temperature≥150◦C sufficed to trigger the viscoplastic buckling of columns(Baˇz ant and Le2008).This conclusion is further supported by simple calculations showing that if,for instance,the column load is raised at temperature250◦C from0.3P t to0.9P t(where P t=failure load=tangent modulus load),the critical time of creep buckling(Baˇz ant and Cedolin2003,chapters8and9)gets shortened from2400hours to1hour(note that,in structural mechanics,the term‘creep buckling’or ‘viscoplastic buckling’represents any time-dependent buckling;on the other hand,in materials science,the term‘creep’is reserved for the time-dependent deformation at stresses<0.5σ0, while the time-dependent deformation at stresses nearσ0is called the‘flow’;Frost and Ashby 1982).Therefore,to decide whether the gravity-driven progressive collapse is the correct explana-tion,the temperature level alone is irrelevant(Baˇz ant and Le2008).It is meaningless and a waste of time to argue about it without calculating the stresses in columns.For low stress, high temperature is necessary to cause collapse,but for high enough stress,even a modestly elevated temperature will cause it.The fact that,after aircraft impact,the loads of some columns must have been close to their strength limit can be clarified by Fig. 1.The asymmetry of aircraft impact damage caused the stiffness centroid of the story to acquire a significant eccentricity,e(Fig.1b).The corresponding bending moment P e of gravity load P=m0g(m0=mass of the initial upper falling part;and g=gravity acceleration)caused nonuniform axial shortening and axial stresses in the surviving columns(Fig.1e,h),which raised the stresses in the columns on the weaker side of story much above the average stress due to gravity.The subsequent heating weakened the overloaded columns on the weaker side(left side in thefigure)more than those on the stronger side,and caused gradually more and more of them to lose their load-carrying capability.This further enlarged e and thus increased the nonuniformity of column deformations and stresses (Fig.1f,i),until the buckling of a sufficient number of columns led to the overall stability loss.The observed multistory inward bowing of some perimeter columns(Fig.1c),which reached 1.40m,must have been a significant factor in stability loss,since the tripling of buckling length reduces the column capacity9-times.As noted by NIST,the bowing must have been caused by sagging of heatedfloor trusses,due to their viscoplastic deformation(however,as pointed out by Baˇz ant and Le2008,the differential thermal expansion could not have been the main cause of bowing since a temperature difference of above1000◦C across thefloor truss would be needed to cause a curvature for which the truss span gets shortened by1.40m).Inabsorbable Kinetic EnergyFirst,let us review the basic argument(Baˇz ant2001;Baˇz ant and Zhou2002).After a drop through at least the height h of one story heated byfire(stage3in Fig.2top),the mass of the upper part of each tower has lost enormous gravitational energy,equal to m0gh.Because the energy dissipation by buckling of the hot columns must have been negligible by comparison, most of this energy must have been converted into kinetic energy K=m0v2/2of the upper part of tower,moving at velocity v.Calculation of energy W c dissipated by the crushing of all columns of the underlying(cold and intact)story showed that,approximately,the kinetic energy of impact K>8.4W c(Eq.3of Baˇz ant and Zhou2002).It is well known that,in inelastic buckling,the deformation must localize into inelastic hinges (Baˇz ant and Cedolin2003,sec.7.10).To obtain an upper bound on W c,the local buckling offlanges and webs,as well as possible steel fracture,was neglected(which means that the ratio K/W c was at least8.4).When the subsequent stories are getting crushed,the loss m0gh of gravitational energy per story exceeds W c exceeds8.4by an ever increasing margin,and so the velocity v of the upper part must increase from one story to the next.This is the basic characteristic of progressive collapse,well known from many previous disasters with causes other thanfire(internal or external explosions,earthquake,lapses in quality control;see,e.g., Levy and Salvadori1992;Baˇz ant and Verdure2007).Merely to get convinced of the inevitability of gravity driven progressive collapse,further analysis is,for a structural engineer,superfluous.Further analysis is nevertheless needed to dispel false myths,and also to acquire full understanding that would allow assessing the danger of progressive collapse in other situations.Generalization of Differential Equation of Progressive CollapseThe gravity-driven progressive collapse of a tower consists of two phases—the crush-down, followed by crush-up(Fig.2bottom),each of which is governed by a different differential equation(Baˇz ant and Verdure2007,pp.312-313).During the crush-down,the falling upper part of tower(C in Fig.2bottom),having a compacted layer of debris at its bottom(zone B), is crushing the lower part(zone A)with negligible damage to itself.During the crush-up,the moving upper part C of tower is being crushed at bottom by the compacted debris B resting on the ground.The fact that the crush-up of entire stories cannot occur simultaneously with the crush-down is demonstrated by the condition of dynamic equilibrium of compacted layer B,along with an estimate of the inertia force of this layer due to vertical deceleration or acceleration;see Eq. 10and Fig.2(f)of Baˇz ant and Verdure(2007).This previous demonstration,however,was only approximate since it did not take into account the variation of crushing forces F c and Fc during the collapse of a story.An accurate analysis of simultaneous(deterministic)crush-upand crush-down is reported in Baˇz ant and Le(2008)and is reviewed in the Appendix,where the differential equations and the initial conditions for a two-way crush are formulated.It is found that,immediately after thefirst critical story collapses,crush fronts will propagate both downwards and upwards.However,the crush-up front will advance into the overlying story only by about1%of its original height h and then stop.Consequently,the effect of the initial two-way crush is imperceptible and the hypothesis that the crush-down and crush-up cannot occur simultaneously is almost exact.The aforementioned distance of initial crush-up would be larger if the column cross sections changed discontinuously right below or right above thefirst collapsed story.However,this does not appear to be the case.A sudden change of column cross section after the crush-down front has advanced by more than a few stories would not produce crush-up because the compacted layer B has already become quite massive and acquired a significant kinetic energy.As an improvement over previous studies,we take into account the fact that a variation of the initial mass density of stories causes the compaction ratioλ(z)to be variable.Consequently, the distribution of velocity v(z)throughout the compacted layer B in the crush-down phase is non-uniform.The equation of motion for the upper part of tower(parts C and B together) during crush-down may be written asd d tz(t)µ(S)˙s(S)d S−gz(t)µ(S)d S=−F c(z,˙z)(1)where t=time,z=vertical(Lagrangian)coordinate=distance of the current crushing front from the initial position of the tower top;the superior dots denote time derivatives;µ(S)= initial specific mass of tower(mass of a story divided by its height)at point of initial coordinate S;˙s(S)=velocity of material point with initial coordinate S.It will suffice to consider the velocity,as well as the momentum density,to be distributed throughout the compacted layer linearly.With these approximations,the crush-down differential equation of motion becomes:d d tm0[1−λ(z)]d zd t+µc l2[2−λ(z)]d zd t−m(z)g=−F c(z,˙z)(crush-down)(2)while the crush-up differential equation of motion has the same form as Eq.17of Baˇz ant and Verdure(2007):m(y)dd t[1−λ(y)]dydt+g=F c(y,˙y)(crush-up)(3)Here Eq.(2)represents a refinement of Eq.12of Baˇz ant and Verdure(2007),while Eq.(3)is identical to their Eq.17because the compacted layer is stationary during crush-up. Furthermore,l=height of compacted layer B,µc=specific mass of compacted layer B per unit height,which is considered to be constant and equal to the maximum possible density of compacted debris;m(z)=cumulative mass of the tower above level z of the crushing front (m(z)=m0+µc l);and F c=resisting force=energy dissipation per unit height;F c(z,˙z)=F b+F s+F a+F e,F b=W d/(1−λ)h(4)where W d(z)=total energy dissipation up to level z,which was assumed by Baˇz ant and Verdure (2007)to consist only of energy F b(per unit height)consumed by buckling of steel columns.In calculations,the largefluctuations of F b as a function of z or y(evident in Figs.3and4 of Baˇz ant and Verdure,2007)are neglected,i.e.,F b is smoothly homogenized.As a refinement of previous analysis,we introduce here a generalization in which we add energy F s(per unit height)consumed by comminution of concretefloor slabs,energy F a required to expel air from the tower,and energy F e required to accelerate the mass of dust and larger fragments ejected from the tower during the impact of upper part;Furthermore,in contrast to previous studies, the compaction ratio will not be assumed as a constant but will be more accurately calculated asλ(z)=(1−κout)µ(z)/µc,andκout=mass shedding fraction=fraction of mass that escapes outside tower perimeter before the end of crush-down(not afterwards).Note that Eq.(2)may be rewritten as[m0(1−λ)+µc l(1−0.5λ)]¨z−mg=−F m−F c,F m=[m0(1−λ)+µc l(1−0.5λ)]˙˙z=¯µ˙z2(5)where F m=force required to accelerate to velocity˙z the stationary mass accreting at the crushing front,and¯µ=d[m0(1−λ)+µc l(1−0.5λ)]/d z=part of the impacted mass per unit height that remains within the tower perimeter.This force causes a greater difference from free fall than do forces F b,F s,F a and F e combined.Upon setting v=˙z,Eq.(2)or(5)was converted to a system of two nonlinearfirst-order differential equations for unknowns v(t)and z(t),which were then integrated numerically with high accuracy using the Runge-Kutta algorithm(note that,for the idealized special case of λ=F c=κout=0and constantµ=d m/d z,Eq.(2)reduces to the differential equation (z˙z)˙=gz,which was formulated and solved byfinite differences by Kausel,2001).As the initial conditions,it is considered that the crushing front initiates at the96th story in the North Tower,and at the81st story in the South Tower(NIST2005).Based on video and seismic records,κout≈0.2,and based on the known typical density of rubble,µc=4.10×106kg/m.These are the values used to calculate the diagrams presented here.However,calculations have also been made for the full range of uncertainty ofκout,as discussed later,and forµc ranging from3.69×106kg/m to4.51×106kg/m.Within these ranges,the changes in the results were not major.Variation of Mass and Buckling Resistance along the HeightNear the top,the specific mass (mass per unit height)µ=1.02×106kg/m.In view of proportionality to the cross section area of columns,µ=1.05×106kg/m at the impact level (81st floor)of South Tower.Although precise data on µ(z )are unavailable,it appears sufficient to use the approximation µ(z )=k 0e k 2z +k 1(where k 0,k 1,k 2=constants),with a smooth transition at the 81st floor to a linear variation all the way down.The condition that H 0µ(z )d z be equal to the total mass of tower (known to be roughly 500,000tons)gives µ=1.46×106kg/m at the base.There are various local complexities whose possible effects were estimated in calculations (e.g.,the fact that 16of 47core columns at the bottom were much more massive than the rest).However,they appeared to have no appreciable effect on the overall response,particularly on the diagram of z (t )and the collapse duration.The total energy dissipation per unit height,which represents the resisting force F c ,consists not only of energy F b dissipated by the inelastic hinges formed during column buckling,but also of energy F s required for comminuting concrete floor slabs,energy F a required for expelling air from the tower,and energy F e required for ejecting particles and fragments.Based on Fig.5and Eq.8of Baˇz ant and Zhou (2002),on Fig.3of Baˇz ant and Verdure (2007),and on Baˇz ant and Cedolin (2003,Sec 8.6),we have,for three-hinge column buckling:F b =u f 0F (u )d u h ,F (u )=N i =12[M a i (θi )+M b i (θ i )]L i sin θi ,θi =arccos 1−u L i (6)where F =axial force resultant of all the columns in the story;u =vertical relative displacement between column ends,u f =final u -value;θ i =2θi ;θi ,θ i =hinge rotations at the ends and middle of column i ,which are functions of u ;M a i ,M b i =bending moments in inelastic hinges at the ends and middle of column i ,as functions of θi or θ i ;and L i =initial clear length of columns i .For plate-type four-hinge buckling (Fig.2c in Baˇz ant and Zhou 2002),similar simple expressions apply.Although some core columns were rectangular,their plastic bending moments M p were nearly proportional to the column cross section areas because,in the weak buckling direction,most core columns had the same width as the perimeter columns.Thus the curve F (u )corre-sponding to perfect plasticity (M a i =M b i =M p i )is not difficult to estimate from the weight of all the columns in a story.However,three effects doubtless intervened to reduce F (u ):1)multi-story buckling of some columns;2)softening due to local plastic flange buckling,and 3)fracture of steel in inelastic hinges (the last two likely occurred only at large buckling deflections for which F (u )is small).The available data are insufficient to make an accurate estimate of these effects,and even the data on the flange thicknesses in the perimeter and core columns of all the stories are missing.So we simply apply to F b an empirical correction factor β(≤1)which is reasonably expected to lie within the range (0.5,0.8)for normal structural steel (yield limit 250MPa),but in the range (0.1,0.3)for the high-strength steel (yield limit 690MPa)which was used for perimeter columns in the lower stories.The high-strength steel has a much lower ductility,which must have caused fractures (with a drop of axial force to zero)very early during buckling,and must have been the cause of formation of large multistory fragments seen to fall from the lower part of tower.Consequently,the energy dissipated (which is equal to the area under the load-displacement curve of column)was probably about the same for high-and normal-strength columns.The effect of uncertainty in F b ,i.e.β=[0.5,0.8],on the collapse time is shown in Fig.8.The remaining diagrams are calculated with β≈2/3for the stories with normal-strength column,while for high-strength columns,βis taken as 0.24,which gives about the same F b as for normal-strength columns.But,for the stories with high-strength perimeter columns,the strength of core columns is 290MPa (i.e.,only slightly higher than for structural steel of normal strength).So we use β≈0.2for these stories and expect the uncertainty range (0.1,0.3).The full range of βwill later be considered in estimating the uncertainty of response.To estimate F b ,we consider that each perimeter column had a box cross section,almost square,with the width of 356mm and web thickness varying from 7.5mm near the top of tower to 125mm near the bottom,while at the aircraft impact level the web was 10mm thick (Wierzbicki and Teng 2003;Hart et al.1985).The core columns had a rectangular box cross section from the tower base up to almost the top,where they had an I-section.The size of corecolumns was not uniform even for a single story.Near the base,some reached1371mm×559 mm,with wall thickness127mm(Engrg.News Record1970).The cross section dimensions can be estimated,albeit only crudely,from recently released architectural blueprints.Based on the area under the buckling curve in Fig.3of Baˇz ant and Verdure(2007),the energy dissipation due to column buckling at the impact zone of the North Tower(96th story) is F b(1−λ)h=0.51GJ(or approximately0.5GJ,as estimated by Baˇz ant and Zhou2002). For other stories,this quantity is scaled according to the approximate cross section area of columns.For the sake of computational effectiveness,and because Eq.(2)is the continuum approximation of a difference equation,the variation of W d/h is considered to be smooth, following a similar equation as that forµ(z).Velocity of Air Ejected from the TowerAn upper bound on area through which the air initially contained within every story gets expelled(Fig.3a)is A w=4ψah c,where4ah c=area of one perimeter wall,a=64m=width of the side of square cross section of tower,h c=3.69m=clear height of one story=distance from the bottom of a story slab to the top of the underlying slab,andψ=vent ratio=ratio of unobstructed(open)area of the perimeter walls to their total area(ψ≤1).The initial mass of air within one story is m a=ρa a2h c,whereρa=1.225kg/m3=mass density of air at atmospheric pressure and room temperature.Just outside the tower perimeter,the air jetting out(Fig.3a)must regain the atmospheric pressure as soon as it exits(White1999,p.149),and its temperature must be roughly equal to the initial temperature(this is a well-known general feature of exhausts,e.g.,from jet engines(White1999,p.149)or pipes(Munson et al.2006)). So,the mass density of exiting airρ≈ρa.The time during which the top slab collapses onto the lower slab≈∆t=h c/˙z=time during which the air is expelled out(which is only about0.07s for stories near the ground). Conservation of the mass of air during the collapse of one story requires thatρA w(v a∆t)=ρV a. Solving this equation gives the average velocity of escaping air just outside the tower perimeter:v a=V aψA w∆t=a˙z4ψh c(7)Since the velocity of the crushing front near the end of North Tower crush-down is,according to the solution of Eq.(2),˙z=47.34m/s(106mph),the velocity of escaping air near the end of crush-down isv a=64m×47.34m/s4ψ×3.69m=205m/s(459mph or0.60Mach)forψ=1340m/s(761mph or1.00Mach)forψ=0.604(8)The vent ratioψ(which is<1)is hard to estimate.It surely varies from story to story,and also during the crushing of one story.Its effective,or average,value could be much less than 1(because some of the perimeter area is doubtless still obstructed early in the crushing of one story,and because much of the air escapes only after the story height has been reduced greatly). In spite of these uncertainties,it is clear that the exit air speed is of the order of500mph and that itsfluctuations must reach the speed of sound.This must,of course,create sonic booms, which are easily mistaken for explosions(supersonic speeds are virtually impossible since the venting would require an orifice shaped similarly to a convergent-divergent nozzle).There are other phenomena that can cause v a to differ from the estimate in Eq.(8).The air pressure surely exhausts the load capacity of thefloor slab for a few microseconds before it is impacted by the layer of compacted debris.So,thefloor slab must crack before the story height is reduced toλh,and the air must begin to leak through the crackedfloor slab into the underlying story,thus increasing the air mass in that story.Obviously some air must also leak into the ceiling which behaves as a porous layer of compacted gravel(it is impossible for the ceiling and thefloor to remainflat and leak no air since otherwise the air pressure would tend to infinity as the ceiling impacts thefloor).All these complex inter-story interactions must cause rapid and large randomfluctuations of internal air pressure and exiting air velocity.On the average,however,what matters is the simple fact that the air must,in one way or another,get expelled from each story of the tower within a very short time interval,which is only0.07s near the end of crush-down of NorthTower.This fact inevitably leads to the average exit velocity estimate in Eq.(7).The high velocity of air jetting out also explains why a large amount of pulverized concrete,drywalls and glass was ejected to a distance of several hundred meters from the tower(Fig.3a).Resisting Forces Due to Ejecting Air and SolidsThe air mass within the confines of one story,which isρa a2h c,gets accelerated from0to velocity v a as it exits the tower perimeter.The kinetic energy acquired by the escaping air of one story just outside the tower perimeter is K a=1ρv a2(a2h c)where a2h c=initial volume of air within the story.The energy dissipated by viscosity offlowing air and by boundary friction is estimated to be negligible.Therefore,virtually all of the kinetic energy of escaping air must be supplied by gravity,and since the spatial derivative of energy is a force(called the material force or configurational force),the vertical resisting force caused by air ejection isF a(z,˙z)=Kh=ρa a432ψ2h c h˙z2(9)Solution of Eq.(2)shows that,at the end of North Tower crush-down,F a≈103MN,which represents about8.01%of the total resisting force F c at the end of crush-down.When thefirst story under the aircraft impact zone gets smashed,F a contributes only2.75%of F c.During the crush-down,the ratio F a/F c starts increasing and reaches the maximum of12.79%when the 83rd story gets impacted,and then decreases due to the increase of F b at lower stories.F a< 5%of F c up to the3rd story crushed,and<10%up to the7th story crushed(by contrast,in building demolitions,which are conducted in the crush-up mode,the crush-up motion begins with zero velocity,and F a<5%of F c for buildings up to about20stories tall,in which case the air resistance can be neglected).The maximum of F a/F c is30.27%,which occurs during the crush-up phase in which F b is small and the velocity high.The average over-pressure of air within the tower is∆p a=F a/a2above the atmospheric pressure,which gives,for the North Tower,7.55kPa(0.075atm),14.31kPa(0.141atm)and 25.15kPa(0.248atm),respectively,during the crushing of the80th story,50th story,and at the end of crush-down,respectively.The last pressure value is enough to break up thefloor slab.The pressure peaks near the end of squeezing of a story are doubtless much higher,as already mentioned,and thus must contribute to the break up of manyfloor slabs(theoretically, the pressure in a thin layer of viscous gas between two colliding parallelflat slabs approaches infinity at the end).The mass that is shed from the tower,characterized byκout,exits at various velocities ranging from nearly0to almost either the air ejection velocity,forfine dust,or to roughly ˙z,for large steel pieces.Instead of complicating our model by some distribution of these velocities,we will simply assume that a certain fraction,κeκout,gets ejected in any direction (horizontal,inclined downward or upward,or almost vertical)at velocity˙z,while the remaining mass(1−κe)κout is shed at nearly vanishing velocity.For a certain empirical value ofκe,this must be energetically equivalent to considering the actual distribution of velocities of ejected solids.As the crushing front advances d z,the mass of solids(dust plus large fragments)that is ejected at velocity˙z isκeκoutµ(z)d z and has kinetic energyκeκoutµ(z)d z(˙z2/2).This must be equal to F e d z,i.e.,to the work of the resisting force F e over distance d z.It follows thatF e=1κeκoutµ(z)˙z2(10) The computation results shown infigures have been run forκe=0.2;however,a broad range ofκe has been considered in computations,as discussed later.For the crush-up,κe must be ignored because the compacted layer is stationary.Energy Dissipated by Comminution(or Fragmentation and Pulver-ization)The energy per unit height,F s,that is dissipated by comminution of concretefloor slabs, can be calculated from the comminution theory,which is a subject well understood by now。

3 2 1
Abstract The main part of this review is devoted to the comprehensive description of coherent radiation by nuclear spins. The theory of nuclear spin superradiance is developed and the experimental observations of this phenomenon are considered. The intriguing problem of how coherence develops from initially incoherent quantum fluctuations is analysed. All main types of coherent radiation by nuclear spins are discussed, which are: free nuclear induction, collective induction, maser generation, pure superradiance, triggered superradiance, pulsing superradiance, punctuated superradiance, and induced emission. The influence of electron-nuclear hyperfine interactions and the role of magnetic anisotropy are studied. Conditions for realizing spin superradiance by magnetic molecules are investigated. The possibility of nuclear matter lasing, accompanied by pion or dibaryon radiation, is briefly touched.

Abstract: A class of fractional order neutral differential equations is considered. By using the method of steps and the theory of differential inequalities, the existence and uniqueness theorems and the finite time stability theorem of the fractional order neutral differential equations are obtained, which extends some corresponding results in [8, 9]. Keywords: existence and uniqueness; finite time stability; fractional order neutral 2010 MR Subject Classification: 34K20; 34K37; 34K40 Document code: A Article ID: 0255-7797(2014)01-0043-08
Vol. 34 ( 2014 ) No. 1
数 学 杂 志
J. of Math. (PRC)
FINITE TIME STABILITY OF FRACTIONAL ORDER NEUTRAL DIFFERENTIAL EQUATIONS
LIU Ke-wei1,2 , JIANG Wei1
(1.School of Mathematical Science, Anhui University, Hefei 230039, China ) (2.School of Mathematics Science, Hefei University of Technology, Hefei 230009, China )