Adiabatic stabilization of excitons in intense terahertz laser

Bull Earthquake Eng(2008)6:645–675DOI10.1007/s10518-008-9078-1ORIGINAL RESEARCH PAPERNumerical analyses of fault–foundation interactionI.Anastasopoulos·A.Callerio·M.F.Bransby·M.C.R.Davies·A.El Nahas·E.Faccioli·G.Gazetas·A.Masella·R.Paolucci·A.Pecker·E.RossignolReceived:22October2007/Accepted:14July2008/Published online:17September2008©Springer Science+Business Media B.V.2008Abstract Field evidence from recent earthquakes has shown that structures can be designed to survive major surface dislocations.This paper:(i)Describes three differentfinite element(FE)methods of analysis,that were developed to simulate dip slip fault rupture propagation through soil and its interaction with foundation–structure systems;(ii)Validates the developed FE methodologies against centrifuge model tests that were conducted at the University of Dundee,Scotland;and(iii)Utilises one of these analysis methods to conduct a short parametric study on the interaction of idealised2-and5-story residential structures lying on slab foundations subjected to normal fault rupture.The comparison between nume-rical and centrifuge model test results shows that reliable predictions can be achieved with reasonably sophisticated constitutive soil models that take account of soil softening after failure.A prerequisite is an adequately refined FE mesh,combined with interface elements with tension cut-off between the soil and the structure.The results of the parametric study reveal that the increase of the surcharge load q of the structure leads to larger fault rupture diversion and“smoothing”of the settlement profile,allowing reduction of its stressing.Soil compliance is shown to be beneficial to the stressing of a structure.For a given soil depthH and imposed dislocation h,the rotation θof the structure is shown to be a function of:I.Anastasopoulos(B)·G.GazetasNational Technical University,Athens,Greecee-mail:ianast@civil.ntua.grA.Callerio·E.Faccioli·A.Masella·R.PaolucciStudio Geotecnico Italiano,Milan,ItalyM.F.BransbyUniversity of Auckland,Auckland,New ZealandM.C.R.Davies·A.El NahasUniversity of Dundee,Dundee,UKA.Pecker·E.RossignolGeodynamique et Structure,Paris,France123(a)its location relative to the fault rupture;(b)the surcharge load q;and(c)soil compliance.Keywords Fault rupture propagation·Soil–structure-interaction·Centrifuge model tests·Strip foundation1IntroductionNumerous cases of devastating effects of earthquake surface fault rupture on structures were observed in the1999earthquakes of Kocaeli,Düzce,and Chi-Chi.However,examples of satisfactory,even spectacular,performance of a variety of structures also emerged(Youd et al.2000;Erdik2001;Bray2001;Ural2001;Ulusay et al.2002;Pamuk et al.2005).In some cases the foundation and structure were quite strong and thus either forced the rupture to deviate or withstood the tectonic movements with some rigid-body rotation and translation but without damage(Anastasopoulos and Gazetas2007a,b;Faccioli et al.2008).In other cases structures were quite ductile and deformed without failing.Thus,the idea(Duncan and Lefebvre1973;Niccum et al.1976;Youd1989;Berill1983)that a structure can be designed to survive with minimal damage a surface fault rupture re-emerged.The work presented herein was motivated by the need to develop quantitative understan-ding of the interaction between a rupturing dip-slip(normal or reverse)fault and a variety of foundation types.In the framework of the QUAKER research project,an integrated approach was employed,comprising three interrelated steps:•Field studies(Anastasopoulos and Gazetas2007a;Faccioli et al.2008)of documented case histories motivated our investigation and offered material for calibration of the theoretical methods and analyses,•Carefully controlled geotechnical centrifuge model tests(Bransby et al.2008a,b)hel-ped in developing an improved understanding of mechanisms and in acquiring a reliable experimental data base for validating the theoretical simulations,and•Analytical numerical methods calibrated against the abovefield and experimental data offered additional insight into the nature of the interaction,and were used in developing parametric results and design aids.This paper summarises the methods and the results of the third step.More specifically: (i)Three differentfinite element(FE)analysis methods are presented and calibratedthrough available soil data.(ii)The three FE analysis methods are validated against four centrifuge experiments con-ducted at the University of Dundee,Scotland.Two experiments are used as a benchmark for the“free-field”part of the problem,and two more for the interaction of the outcrop-ping dislocation with rigid strip foundations.(iii)One of these analysis methods is utilised in conducting a short parametric study on the interaction of typical residential structures with a normal fault rupture.The problem studied in this paper is portrayed in Fig.1.It refers to a uniform cohesionless soil deposit of thickness H at the base of which a dip-slip fault,dipping at angle a(measured from the horizontal),produces downward or upward displacement,of vertical component h.The offset(i.e.,the differential displacement)is applied to the right part of the model quasi-statically in small consecutive steps.123hx O:“f o c u s ”O ’:“e p i c e n t e r ”Hanging wallFootwallyLW –LW hx O:“fo c u s ”O ’:“e p i c e n t e r ”Hanging wallFootwallyL W –LWq BStrip Foundation s(a )(b)Fig.1Definition and geometry of the studied problem:(a )Propagation of the fault rupture in the free field,and (b )Interaction with strip foundation of width B subjected to uniform load q .The left edge of the foundation is at distance s from the free-field fault outcrop2Centrifuge model testingA series of centrifuge model tests have been conducted in the beam centrifuge of the University of Dundee (Fig.2a)to investigate fault rupture propagation through sand and its in-teraction with strip footings (Bransby et al.2008a ,b ).The tests modelled soil deposits of depth H ranging from 15to 25m.They were conducted at accelerations ranging from 50to 115g.A special apparatus was developed in the University of Dundee to simulate normal and reverse faulting.A central guidance system and three aluminum wedges were installed to impose displacement at the desired dip angle.Two hydraulic actuators were used to push on the side of a split shear box (Fig.2a)up or down,simulating reverse or normal faulting,respectively.The apparatus was installed in one of the University of Dundee’s centrifuge strongboxes (Fig.2b).The strongbox contains a front and a back transparent Perspex plate,through which the models are monitored in flight.More details on the experimental setup can be found in Bransby et al.(2008a ).Displacements (vertical and horizontal)at different123Fig.2(a)The geotechnicalcentrifuge of the University ofDundee;(b)the apparatus for theexperimental simulation of faultrupture propagation through sandpositions within the soil specimen were computed through the analysis of a series of digital images captured as faulting progressed using the Geo-PIV software(White et al.2003).Soil specimens were prepared within the split box apparatus by pluviating dry Fontainebleau sand from a specific height with controllable massflow rate.Dry sand samples were prepared at relative densities of60%.Fontainebleau sand was used so that previously published laboratory element test data(e.g Gaudin2002)could be used to select drained soil parameters for thefinite element analyses.The experimental simulation was conducted in two steps.First,fault rupture propagation though soil was modelled in the absence of a structure(Fig.1a),representing the free-field part of the problem.Then,strip foundations were placed at a pre-specified distance s from the free-field fault outcrop(Fig.1b),and new tests were conducted to simulate the interaction of the fault rupture with strip foundations.3Methods of numerical analysisThree different numerical analysis approaches were developed,calibrated,and tested.Three different numerical codes were used,in combination with soil constitutive models ranging from simplified to more sophisticated.This way,three methods were developed,each one corresponding to a different level of sophistication:(a)Method1,using the commercial FE code PLAXIS(2006),in combination with a simplenon-associated elastic-perfectly plastic Mohr-Coulomb constitutive model for soil; 123Foundation : 2-D Elastic Solid Elements Elastic BeamElementsInterfaceElements hFig.3Method 1(Plaxis)finite element diecretisation(b)Method 2,utilising the commercial FE code ABAQUS (2004),combined with a modifiedMohr-Coulomb constitutive soil model taking account of strain softening;and(c)Method 3,making use of the FE code DYNAFLOW (Prevost 1981),along with thesophisticated multi-yield constitutive model of Prevost (1989,1993).Centrifuge model tests that were conducted in the University of Dundee were used to validate the effectiveness of the three different numerical methodologies.The main features,the soil constitutive models,and the calibration procedure for each one of the three analysis methodologies are discussed in the following sections.3.1Method 13.1.1Finite element modeling approachThe first method uses PLAXIS (2006),a commercial geotechnical FE code,capable of 2D plane strain,plane stress,or axisymmetric analyses.As shown in Fig.3,the finite element mesh consists of 6-node triangular plane strain elements.The characteristic length of the elements was reduced below the footing and in the region where the fault rapture is expected to propagate.Since a remeshing technique (probably the best approach when dealing with large deformation problems)is not available in PLAXIS ,at the base of the model and near the fault starting point,larger elements were introduced to avoid numerical inaccuracies and instability caused by ill conditioning of the element geometry during the displacement application (i.e.node overlapping and element distortion).The foundation system was modeled using a two-layer compound system,consisting of (see Fig.3):•The footing itself,discretised by very stiff 2D elements with linear elastic behaviour.The pressure applied by the overlying building structure has been imposed to the models through the self weight of the foundation elements.123Fig.4Method1:Calibration of constitutive model parameters utilising the FE code Tochnog;(a)oedometer test;(b)Triaxial test,p=90kPa•Beam elements attached to the nodes at the bottom of the foundation,with stiffness para-meters lower than those of the footing to avoid a major stiffness discontinuity between the underlying soil and the foundation structure.•The beam elements are connected to soil elements through an interface with a purely frictional behaviour and the same friction angleϕwith the soil.The interface has a tension cut-off,which causes a gap to develop between soil and foundation in case of detachment. Due to the large imposed displacement reached during the centrifuge tests(more than3m in several cases),with a relative displacement of the order of10%of the modeled soil height, the large displacement Lagrangian description was adopted.After an initial phase in which the geostatic stresses were allowed to develop,the fault displacement has been monotonically imposed both on the right side and the right bottom boundaries,while the remaining boundaries of the model have beenfixed in the direction perpendicular to the side(Fig.3),so as to reproduce the centrifuge test boundary conditions.3.1.2Soil constitutive model and calibrationThe constitutive model adopted for all of the analyses is the standard Mohr-Coulomb for-mulation implemented in PLAXIS.The calibration of the elastic and strength parameters of the soil had been conducted during the earlier phases of the project by means of the FEM code Tochnog(see the developer’s home page ),adopting a rather refined and user-defined constitutive model for sand.This model was calibrated with a set of experimental data available on Fontainebleau sand(Gaudin2002).Oedometer tests (Fig.4a)and drained triaxial compression tests(Fig.4b)have been simulated,and sand model parameters were calibrated to reproduce the experimental results.The user-defined model implemented in Tochnog included a yielding function at the critical state,which corresponds to the Mohr-Coulomb failure criterion.A subset of those parameters was then utilised in the analysis conducted using the simpler Mohr-Coulomb model of PLAXIS:•Angle of frictionϕ=37◦•Young’s Modulus E=675MPa•Poisson’s ratioν=0.35•Angle of Dilationψ=0◦123hFoundation : Elastic Beam ElementsGap Elements Fig.5Method 2(Abaqus)finite element diecretisationThe assumption of ψ=0and ν=0.35,although not intuitively reasonable,was proven to provide the best fit to experimental data,both for normal and reverse faulting.3.2Method 23.2.1Finite element modeling approachThe FE mesh used for the analyses is depicted in Fig.5(for the reverse fault case).The soil is now modelled with quadrilateral plane strain elements of width d FE =1m.The foun-dation,of width B ,is modelled with beam elements.It is placed on top of the soil model and connected through special contact (gap)elements.Such elements are infinitely stiff in compression,but offer no resistance in tension.In shear,their behaviour follows Coulomb’s friction law.3.2.2Soil constitutive modelEarlier studies have shown that soil behaviour after failure plays a major role in problems related to shear-band formation (Bray 1990;Bray et al.1994a ,b ).Relatively simple elasto-plastic constitutive models,with Mohr-Coulomb failure criterion,in combination with strain softening have been shown to be effective in the simulation of fault rupture propagation through soil (Roth et al.1981,1982;Loukidis 1999;Erickson et al.2001),as well as for modelling the failure of embankments and slopes (Potts et al.1990,1997).In this study,we apply a similar elastoplastic constitutive model with Mohr-Coulomb failure criterion and isotropic strain softening (Anastasopoulos 2005).Softening is introduced by reducing the mobilised friction angle ϕmob and the mobilised dilation angle ψmob with the increase of plastic octahedral shear strain:123ϕmob=ϕp−ϕp−ϕresγP fγP oct,for0≤γP oct<γP fϕres,forγP oct≥γP f(1)ψmob=⎧⎨⎩ψp1−γP octγP f,for0≤γP oct<γP fψres,forγP oct≥γP f⎫⎬⎭(2)whereϕp andϕres the ultimate mobilised friction angle and its residual value;ψp the ultimate dilation angle;γP f the plastic octahedral shear strain at the end of softening.3.2.3Constitutive model calibrationConstitutive model parameters are calibrated through the results of direct shear tests.Soil response can be divided in four characteristic phases(Anastasopoulos et al.2007):(a)Quasi-elastic behavior:The soil deforms quasi-elastically(Jewell and Roth1987),upto a horizontal displacementδx y.(b)Plastic behavior:The soil enters the plastic region and dilates,reaching peak conditionsat horizontal displacementδx p.(c)Softening behavior:Right after the peak,a single horizontal shear band develops(Jewelland Roth1987;Gerolymos et al.2007).(d)Residual behavior:Softening is completed at horizontal displacementδx f(δy/δx≈0).Then,deformation is accumulated along the developed shear band.Quasi-elastic behaviour is modelled as linear elastic,with secant modulus G S linearly incre-asing with depth:G S=τyγy(3)whereτy andγy:the shear stress and strain atfirst yield,directly measured from test data.After peak conditions are reached,it is assumed that plastic shear deformation takes placewithin the shear band,while the rest of the specimen remains elastic(Shibuya et al.1997).Scale effects have been shown to play a major role in shear localisation problems(Stone andMuir Wood1992;Muir Wood and Stone1994;Muir Wood2002).Given the unavoidableshortcomings of the FE method,an approximate simplified scaling method(Anastasopouloset al.2007)is employed.The constitutive model was encoded in the FE code ABAQUS(2004).Its capability toreproduce soil behaviour has been validated through a series of FE simulations of the directshear test(Anastasopoulos2005).Figure6depicts the results of such a simulation of denseFontainebleau sand(D r≈80%),and its comparison with experimental data by Gaudin (2002).Despite its simplicity and(perhaps)lack of generality,the employed constitutivemodel captures the predominant mode of deformation of the problem studied herein,provi-ding a reasonable simplification of complex soil behaviour.3.3Method33.3.1Finite element modeling approachThefinite element model used for the analyses is shown for the normal fault case in Fig.7.The soil is modeled with square,quadrilateral,plane strain elements,of width d FE=0.5m. 123Fig.6Method 2:Calibration ofconstitutive model—comparisonbetween laboratory direct sheartests on Fontainebleau sand(Gaudin 2002)and the results ofthe constitutive modelx D v3.3.2Soil constitutive ModelThe constitutive model is the multi-yield constitutive model developed by Prevost (1989,1993).It is a kinematic hardening model,based on a relatively simple plasticity theory (Prevost 1985)and is applicable to both cohesive and cohesionless soils.The concept of a “field of work-hardening moduli”(Iwan 1967;Mróz 1967;Prevost 1977),is used by defining a collection f 0,f 1,...,f n of nested yield surfaces in the stress space.V on Mises type surfaces are employed for cohesive materials,and Drucker-Prager/Mohr-Coulomb type surfaces are employed for frictional materials (sands).The yield surfaces define regions of constant shear moduli in the stress space,and in this manner the model discretises the smooth elastic-plastic stress–strain curve into n linear segments.The outermost surface f n represents a failure surface.In addition,accounting for experimental evidence from tests on frictional materials (de 1987),a non-associative plastic flow rule is used for the dilatational component of the plastic potential.Finally,the material hysteretic behavior and shear stress-induced anisotropic effects are simulated by a kinematic rule .Upon contact,the yield surfaces are translated in the stress space by the stress point,and the direction of translation is selected such that the yield surfaces do not overlap,but remain tangent to each other at the stress point.3.3.3Constitutive model parametersThe required constitutive parameters of the multi-yield constitutive soil model are summari-sed as follows (Popescu and Prevost 1995):a.Initial state parameters :mass density of the solid phase ρs ,and for the case of porous saturated media,porosity n w and permeability k .b.Low strain elastic parameters :low strain moduli G 0and B 0.The dependence of the moduli on the mean effective normal stress p ,is assumed to be of the following form:G =G 0 p p 0 n B =B 0 p p 0n (4)and is accounted for,by introducing two more parameters:the power exponent n and the reference effective mean normal stress p 0.c.Yield and failure parameters :these parameters describe the position a i ,size M i and plastic modulus H i ,corresponding to each yield surface f i ,i =0,1,...n .For the case of pressure sensitive materials,a modified hyperbolic expression proposed by Prevost (1989)and Griffiths and Prévost (1990)is used to simulate soil stress–strain relations.The necessary parameters are:(i)the initial gradient,given by the small strain shear modulus G 0,and (ii)the stress (function of the friction angle at failure ϕand the stress path)and strain,εmax de v ,levels at failure.Hayashi et al.(1992)improved the modified hyperbolic model by introducing a new parameter—a —depending on the maximum grain size D max and uniformity coefficient C u .Finally,the coefficient of lateral stress K 0is necessary to evaluate the initial positions a i of the yield surfaces.d.Dilation parameters :these are used to evaluate the volumetric part of the plastic potentialand consist of:(i)the dilation (or phase transformation)angle ¯ϕ,and (ii)the dilation parameter X pp ,which is the scale parameter for the plastic dilation,and depends basically on relative density and sand type (fabric,grain size).With the exception of the dilation parameter,all the required constitutive model parameters are traditional soil properties,and can be derived from the results of conventional laboratory 123Table1Constitutive model parameters used in method3Number of yield surfaces20Power exponent n0.5Shear modulus G at stress p1 (kPa)75,000Bulk modulus at stress p1(kPa)200,000Unit massρ(t.m−3) 1.63Cohesion0 Reference mean normal stressp1(kPa)100Lateral stress coefficient(K0)0.5Dilation angle in compression (◦)31Dilation angle in extension(◦)31Ultimate friction angle in compression(◦)41.8Ultimate friction angle inextension(◦)41.8Dilation parameter X pp 1.65Max shear strain incompression0.08Max shear strain in extension0.08Generation coefficient in compressionαc 0.098Generation coefficient inextensionαe0.095Generation coefficient in compressionαlc 0.66Generation coefficient inextensionαle0.66Generation coefficient in compressionαuc 1.16Generation coefficient inextensionαue1.16(e.g.triaxial,simple shear)and in situ(e.g.cone penetration,standard penetration,wave velocity)soil tests.The dilational parameter can be evaluated on the basis of results of liquefaction strength analysis,when available;further details can be found in Popescu and Prevost(1995)and Popescu(1995).Since in the present study the sand material is dry,the cohesionless material was modeled as a one-phase material.Therefore neither the soil porosity,n w,nor the permeability,k,are needed.For the shear stress–strain curve generation,given the maximum shear modulus G1,the maximum shear stressτmax and the maximum shear strainγmax,the following functional relationship has been chosen:For y=τ/τmax and x=γ/γr,withγr=τmax/G1,then:y=exp(−ax)f(x,x l)+(1−exp(−ax))f(x,x u)where:f(x,x i)=(2x/x i+1)x i−1/(2x/x i+1)x i+1(5)where a,x l and x u are material parameters.For further details,the reader is referred to Hayashi et al.(1992).The constitutive model is implemented in the computer code DYNAFLOW(Prevost1981) that has been used for the numerical analyses.3.3.4Calibration of model constitutive parametersTo calibrate the values of the constitutive parameters,numerical triaxial tests were simulated with DYNAFLOW at three different confining pressures(30,60,90kPa)and compared with the results of available physical tests conducted on the same material at the same confining pressures.The parameters are defined based on the shear stress versus axial strain curve and volumetric strain versus axial strain curve.Figure8illustrates the comparisons between numerical simulations and physical tests in terms of volumetric strain and shear stress versus123Table2Summary of main attributes of the centrifuge model testsTest Faulting B(m)q(kPa)s(m)g-Level a D r(%)H(m)L(m)W(m)h max(m) 12Normal Free—field11560.224.775.723.53.1528Reverse Free—field11560.815.175.723.52.5914Normal10912.911562.524.675.723.52.4929Reverse10919.211564.115.175.723.53.30a Centrifugal accelerationFig.9Test12—Free-field faultD r=60%Fontainebleau sand(α=60◦):Comparison ofnumerical with experimentalvertical displacement of thesurface for bedrock dislocationh=3.0m(Method1)and2.5m(Method2)[all displacements aregiven in prototype scale]Structure Interaction(FR-SFSI):(i)Test14,normal faulting at60◦;and(ii)Test29,reverse faulting at60◦.In this case,the comparison is conducted for all of the developed numerical analysis approaches.The main attributes of the four centrifuge model tests used for the comparisons are syn-opsised in Table2,while more details can be found in Bransby et al.(2008a,b).4.1Free-field fault rupture propagation4.1.1Test12—normal60◦This test was conducted at115g on medium-loose(D r=60%)Fontainebleau sand,simu-lating normal fault rupture propagation through an H=25m soil deposit.The comparison between analytical predictions and experimental data is depicted in Fig.9in terms of vertical displacement y at the ground surface.All displacements are given in prototype scale.While the analytical prediction of Method1is compared with test data for h=3.0m,in the case of Method2the comparison is conducted at slightly lower imposed bedrock displacement: h=2.5m.This is due to the fact that the numerical analysis with Method2was conducted without knowing the test results,and at that time it had been agreed to set the maximum displacement equal to h max=2.5m.However,when test results were publicised,the actually attained maximum displacement was larger,something that was taken into account in the analyses with Method1.As illustrated in Fig.9,Method2predicts almost correctly the location of fault out-cropping,at about—10m from the“epicenter”,with discrepancies limited to1or2m.The deformation can be seen to be slightly more localised in the centrifuge test,but the comparison between analytical and experimental shear zone thickness is quite satisfactory.The vertical displacement profile predicted by Method1is also qualitatively acceptable.However,the123Method 2Centrifuge Model TestR1S1Method 1(a )(b)(c)Fig.10Test 12—-Normal free-field fault rupture propagation through H =25m D r =60%Fontainebleau sand:Comparison of (a )Centrifuge model test image,compared to FE deformed mesh with shear strain contours of Method 1(b ),and Method 2(c ),for h =2.5mlocation of fault rupture emergence is a few meters to the left compared with the experimen-tal:at about 15m from the “epicenter”(instead of about 10m).In addition,the deformation predicted by Method 1at the ground surface computed using method 1is widespread,instead of localised at a narrow band.FE deformed meshes with superimposed shear strain contours are compared with an image from the experiment in Fig.10,for h =2.5m.In the case of Method 2,the comparison can be seen to be quite satisfactory.However,it is noted that the secondary rupture (S 1)that forms in the experiment to the right of the main shear plane (R 1)is not predicted by Method 2.Also,experimental shear strain contours (not shown herein)are a little more diffuse than the FE prediction.Overall,the comparison is quite satisfactory.In the case of Method 1,the quantitative details are not in satisfactory agreement,but the calculation reveals a secondary rupture to the right of the main shear zone,consistent with the experimental image.4.1.2Test 28—reverse 60◦This test was also conducted at 115g and the sand was of practically the same relative density (D r =61%).Given that reverse fault ruptures require larger normalised bedrock123Fig.11Test28—Reversepropagation through H=15mD r=60%Fontainebleau sand:Comparison of numerical withexperimental verticaldisplacement of the surface forbedrock dislocation h=2.0m(all displacements are given inprototype scale)displacement h/H to propagate all the way to the surface(e.g.Cole and Lade1984;Lade et al.1984;Anastasopoulos et al.2007;Bransby et al.2008b),the soil depth was set at H=15m.This way,a larger h/H could be achieved with the same actuator.Figure11compares the vertical displacement y at the ground surface predicted by the numerical analysis to experimental data,for h=2.0m.This time,both models predict correctly the location of fault outcropping(defined as the point where the steepest gradient is observed).In particular,Method1achieves a slightly better prediction of the outcropping location:−10m from the epicentre(i.e.,a difference of1m only,to the other direction). Method2predicts the fault outbreak at about−7m from the“epicenter”,as opposed to about −9m of the centrifuge model test(i.e.,a discrepancy of about2m).Figure12compares FE deformed meshes with superimposed shear strain contours with an image from the experiment,for h=2.5m.In the case of Method2,the numerical analysis seems to predict a distinct fault scarp,with most of the deformation localised within it.In contrast,the localisation in the experiment is clearly more intense,but the fault scarp at the surface is much less pronounced:the deformation is widespread over a larger area.The analysis with Method1is successful in terms of the outcropping location.However,instead of a single rupture,it predicts the development of two main ruptures(R1and R2),accompanied by a third shear plane in between.Although such soil response has also been demonstrated by other researchers(e.g.Loukidis and Bouckovalas2001),in this case the predicted multiple rupture planes are not consistent with experimental results.4.2Interaction with strip footingsHaving validated the effectiveness of the developed numerical analysis methodologies in simulating fault rupture propagation in the free-field,we proceed to the comparisons of experiments with strip foundations:one for normal(Test14),and one for reverse(Test29) faulting.This time,the comparison is extended to all three methods.4.2.1Test14—normal60◦This test is practically the same with the free-field Test12,with the only difference being the presence of a B=10m strip foundation subjected to a bearing pressure q=90kPa.The foundation is positioned so that the free-field fault rupture would emerge at distance s=2.9m from the left edge of the foundation.123。

A study on the mean crushing strength of hexagonalmulti-cell thin-walled structuresZhonghao Bai a,Hourui Guo a,Binhui Jiang a,n,Feng Zhu b,Libo Cao aa The State Key Laboratory of Advanced Design and Manufacturing for Vehicle Body,Hunan University,Hunan410082,Chinab Bioengineering Center,Wayne State University,MI48201,USAa r t i c l e i n f oArticle history:Received6December2013Received in revised form21February2014Accepted23February2014Available online20March2014Keywords:Multi-cell thin-walled structuresAluminum honeycombQuasi-static compressionMean crushing strengtha b s t r a c tThis paper presents a study to extend a previously developed theoretical model to predict the crushingbehavior of hexagonal multi-cell thin-walled structures,e.g.honeycombs under quasi-static loading.Thelow speed compressive tests were conducted on three types of aluminum honeycomb panels.Based onthe test data and existing theoretical models,a new analytical model was developed to predict its meancrushing strength.Some key parameters in this new model were determined with thefinite element(FE)method.Then the predictions based on the new model were compared with the results reported in thepublished literature.It has been shown that the new model has a similar or better performancecompared to its counterparts.Considering its concise expression,the newly developed model can bedeemed as a convenient computational tool in engineering practice.&2014Elsevier Ltd.All rights reserved.1.IntroductionHexagonal multi-cell thin-walled structures have been widely usedas lightweight energy absorbers to improve the crashworthiness andshock resistance of aircrafts,ships and vehicles[1–4].To evaluate theenergy absorption performance of such structures,the mean crushingstrength should be calculatedfirst,as it determines the energydissipated until the material is compacted.Numerous studies havebeen conducted towards developing theoretical models to predict themean crushing strength for thin-walled structures.For instance,McFarland[5]developed a semi-empirical model to predict the meancrushing strength of hexagonal cell structures.Gibson et al.[6]andGibson and Ashby[7]developed another theoretical model forhoneycombs,which was then validated against experimental data.De Oliveira and Wierzbicki[8],Wierzbicki[9],Wierzbicki andAbramowicz[10]introduced a super folding theory for predictingthe mean crushing strength of thin-walled structures.Wu and Jiang[11]compared the predictions based on the theory by Wierzbicki andco-workers against the quasi-static compression test results onhoneycombs.They found that the theoretical predictions under-estimated the experimental data.Chen and Wierzbicki[12]simplifiedthe super folding theory and improved its performance.This simpli-fied model was then adopted by Zhang et al.[13]to successfullypredict the mean crushing strength of multi-cell square columns.ZareiMahmoudabadi and Sadighi[14–16]further improved Wierzbicki'ssuper folding element theory by considering a more detailed geo-metric change during the structural deformation.Although a betteraccuracy was obtained,the equations have many parameters and theyare not convenient to use in practice.In this paper,a typical hexagonal multi-cell thin-walled struc-ture,i.e.honeycombs were tested in the out-of-plane directionunder quasi-static loading conditions.It was found that adhesivefailure has significant effect on the folding behavior.However,such effect on energy dissipation was not considered in existingtheoretical derivations.The simplified super folding elementtheory by Chen and Wierzbicki[12]was extended to derive theequation to calculate the mean crushing strength of hexagonalhoneycombs.Double-thickness wall adhesive failure effect wasconsidered in model mercial nonlinearfiniteelement(FE)code ABAQUS/EXPLICIT(Version6.10,3DS SIMULIA,RI,USA)was employed to simulate quasi-static compressiveresponse of honeycombs with different materials,foil thicknessesand cell sizes.Based on the simulation results,several parametersin the new formula for hexagonal aluminum honeycomb weredetermined.Finally,the performance of the new formula wascompared with the predictions based on the theories by McFar-land,Wierzbicki and Zarei Mahmoudabadi and Sadighi.2.Experiment studyQuasi-static compression tests were conducted on three typesof aluminum honeycomb samples using an Instron machine.Theconfigurations of the samples were different in materials and cellContents lists available at ScienceDirectjournal homepage:/locate/twsThin-Walled Structures/10.1016/j.tws.2014.02.0240263-8231/&2014Elsevier Ltd.All rightsreserved.n Corresponding author.Tel.:þ8613787054434;fax:þ8673188823715.E-mail address:jjhhzz123@(B.Jiang).Thin-Walled Structures80(2014)38–45sizes,as detailed in Table 1.The dimensions of specimen were L :150mm ÂW :150mm ÂH :50mm (shown in Fig.1).Because the mean crushing strength was the main focus in this study,the specimens were not compressed to the densi fication stage.The typical nominal stress –strain curves and the folding pattern of the three aluminum honeycomb samples are shown in Figs.2and 3,respectively.The mean crushing strength s m of each specimen was obtained by integrating the nominal stress s with respect to nominal strain εand then divided by the integrationstrain range (0.2–0.7),as shown in Eq.(1).Each type sample was tested three times,and the average mean crushing strengths were 0.240MPa,1.184MPa and 0.398MPa.s m ¼1ε2Àε1Z ε2ε1s d εð1Þwhere ε1and ε2are the start and end point nominal strain of the densi fication stage,respectively;and s is the nominal stress of densi fication.During the compression,the folding pattern of the cells at the sample's boundary was irregular due to lack of constraints and thus it was not considered in the subsequent analysis.To inves-tigate the crushing behavior in more detail,a piece of Type 2aluminum honeycomb sample after test was taken and it was flattened for better examination as shown in Fig.4a.Gibson and Ashby [7]pointed out that various deformation modes may be observed during the honeycomb out-of-plane compression,including linear elastic deformation,elastic buckling,plastic deformation and brittle fracture.In the current study,brittle deformation was not seen as the base material,i.e.aluminum is a ductile material.Since aluminum honeycombs were manufac-tured with expended/corrugated and adhesive [17,18],bending of both single-layered/double-layered wall and adhesive failure of double-thickness wall should be considered.In Fig.4a,adhesive failure of double-layered wall and plastic hinge lines were clearly observed.The honeycomb was assumed to be made of an ideal rigid-plastic material [9].The flattened cell walls generally consist of two types of planes:(1)pentagonal planes and triangular planes due to the adhesive failure (Fig.4a).Fig.4b is a sketch showing the detailed folding pattern,where the green transparent area repre-sents an area with double-layered folding;the yellow area is a triangular plane which is dominated by adhesive failure;the khaki and blue areas are pentagonal planes and triangular planes on the single-layered parison of all three type samples indi-cated that the area of yellow triangular plane had a great effect on the width of the adhesive failure region,and the height of yellow triangular plane was equal to the width of the adhesive failure region.Fig.4c shows a paper model of folding process (the region enclosed by dotted line was the flattened area shown in Fig.4a).3.Theory3.1.Super folding element theorySuper folding element theory was proposed by Wierzbicki [9].In this theory,the cell was assumed to be made of an ideal rigid-plastic material [9].Fig.5shows a cell which includes four trapezoidal planes moving as rigid bodies (areas 1–4)and two sections of curved surfaces.Based on the kinematic analysis of a basic folding element during crushing,the work of external force is equal to the energy dissipation by toroid shell extension and plastic hinge lines formation of on the curved surfaces.According to energy conservation law,the mean compressive crushing strength could be obtained using the following equation [9,10]:s m ¼16:56s 0hs 5=3ð2Þwhere s 0is the flow stress of base material;s is the minor cell diameter which is the distance between the two parallel cell edges;and h is the cell wall thickness.Chen and Wierzbicki [12]simpli fied the original super folding element theory in [9].In the new theory,the total energy dissipation is separated into two parts:dissipation during the plastic hinge lines formation and during plastic flow.The simpli-fied super folding element theory can be conveniently used toTable 1Three types of aluminum honeycomb samples tested.Type numberMaterialWallthickness,t (mm)Length of cell edge,D (mm)Nominal density (kg/m 3)Type 1Al30030.0711.0326.3Type 2Al30030.07 4.0473.1Type 3Al50520.0711.0326.3Fig.1.A honeycomb specimen of Type 1for compression tests.0.00.10.20.30.40.50.60.70.80.00.20.40.60.81.01.21.41.61.82.0N o m i n a l S t r e s s (M P a )Nominal StrainFig.2.Typical nominal stress –strain curves of three types of aluminum honeycomb samples under quasi-static compression.Z.Bai et al./Thin-Walled Structures 80(2014)38–4539analyze complicated multi-cell thin-walled structures[12].For instance,Zhang et al.[13]used the simplified super folding element theory to study the multi-cell square thin-walled struc-tures which was divided into three basic components:the corner part,the crisscross part and T-shape part.The mean crushing force P m can be obtained asP m¼s0tffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiðN cþ4N0þ2N TÞπL c tpð3Þwhere N c,N0,and N T denote the number of corner part,crisscrosspart and T-shape part,respectively;and L c t is the crosssectional area.3.2.Theoretical study of mean crushing stress of hexagonalaluminum honeycombA“Y”unit with a length2H(Fig.6)is typically used as the basicfolding element to analyze the out-of-plane behavior of hexagonalhoneycombs in theoretical or FE analysis[8,15,16,19].According toexperimental observation,aluminum foil deformation is mainlygoverned by extension and bending of the planar trapezoidalelements around plastic hinge.The adhesive joining the double-layered wall was separated in part,and thus affecting the defor-mation behavior of basic element and energy dissipation.Tosimplify the analysis,the effect of adhesive failure in energydissipation was not considered in this study[9].According tosimplified super folding element theory,the total energy absorp-tion U of a Y unit after completely folding consists of the energydissipated by plastic hinge lines formation E h and the energydissipated by plasticflow E f.According to Wierzbicki's super folding element theory andChen's simplified theory[9,12],energy absorption duringthe Fig.3.The typical folding patterns of three types of aluminumhoneycomb.Fig.4.Deformation mode of Type2aluminum honeycomb after compression:(a)flattened face of single-thickness wall;(b)a sketch of the deformation pattern;and(c)paper model of folding process.(For interpretation of the references to color in thisfigure,the reader is referred to the web version of thisarticle.)Fig.5.Basic folding element[7].Z.Bai et al./Thin-Walled Structures80(2014)38–4540plastic hinge line formation is dependent on the hinge line lengthand the bending angle of the rigid body around hinge line.In Fig.7(a),dashed line AB formed a double-layered plastic hinge line,andline BC formed two single-layered plastic hinge lines due toadhesive failure in the folding procedure.Therefore the energydissipated E h1on the double-layered hinge line E h1is mainlyinduced by hinge line AB and hinge lines IJ and MN,which areclamped upper and lower boundary.The rotational angle of hingeline AB isθ1¼πwhile the angle for IJ and MN isθ2¼π=2. According to folding pattern observed,the length of the adhesivefailure region is equal to height of the yellow triangular plane(Fig.4b).The length of the triangular plane is proportional to thehalf-wavelength[12].Therefore adhesive failure length L gf is alsoproportional to half-wavelength H.Their proportional coefficientis denoted asα,i.e.L gf¼αH.The length of hinge line AB,IJ and MN L1is equal to CÀαH(where C is the length of Y unit in Fig.6a). Then the energy dissipation by double-layered plastic hinge lines of a Y unit can be determined as follows:E h1¼M0dθ1L1þ2M0dθ2L1ð4Þwhere M0d¼1s0ð2tÞ2being the plastic bending moment;t is the foil thickness;and s0is the effectiveflow stress of the base material which is defined in the following equation[16]:s0¼ffiffiffiffiffiffiffiffiffiffiffis y s u1þnrð5Þwhere s y and s u denote the yield strength and the ultimatestrength of the base material,respectively;and n is the exponentof the power law.The single-layered hinge lines of a Y unit include hinges BC,BEand BG,which are formed by folding in the adhesive failure area aswell as hinge lines CD,EF,GH,LQ and OP on the single-layeredwall(Fig.7(a)).The energy dissipated by hinge lines JB and BN wasnot considered due to the quite small rotation during the crushingprocess.The length of line BC is equal to the length of adhesivefailure area L gf¼αH.The length of line BE,EF,BG and GH L2isproportional to the half-wavelength H[12],i.e.L2¼βH,whereβisthe proportional coefficient.The length of CD,LQ and OP is equalto the side length of Y unit C.The angle of hinge lines LQ and OP isθ2¼π=2,and the angles for other lines are identical,i.e.θ1¼π.Sothe energy absorption by single-layered plastic hinge lines of a Yunit can be obtained by the following equation:E h2¼M0θ1ðCþL gfþ4L2Þþ2M0θ2ðCþL gfÞð6ÞwhereM0¼1=4s0t2:The total energy absorption by plastic hinge lines of a basic Y unitduring the folding process isE h¼E h1þ2E h2¼8M0πðCÀαHÞþ4M0πCþ4αM0πHþ8βM0πH¼12M0πCþð8βÀ4αÞM0πHð7ÞFig.6.Sketch showing the Y-section and a basic folding element in the hexagonal honeycomb:(a)Y-section to be analyzed;and(b)a basic foldingelement.Fig.7.Hinge lines and plasticflowing region:(a)Hinge lines;and(b)plasticflowing region.Z.Bai et al./Thin-Walled Structures80(2014)38–4541where αis the proportional coef ficient between adhesive failure length and half-wavelength H ;and βis the proportional coef fi-cient between plastic hinge line length L 2and half-wavelength H .In the simpli fied super folding element theory,the area of plastic flowing region of a Y unit is equal to H 2/2[12].Considering the uncertainty of the folding process,here we assume the coef ficient between plastic region area S and H 2as γ.The sketch of plastic flowing region is shown in Fig.7(b).In the folding process of a Y unit,there are two single-layered plastic flowing regions with the area γH 2.The energy absorption due to mem-brane effect in one folding element can be evaluated by integrating the membrane stress s 0t in plastic flowing area [9].So the energy dissipated due to plastic flowing can be calculated by the following equation:E f ¼2ZS ðs 0t ÞdS ¼2γs 0tH 2¼8γM 0H 2ð8ÞThen the total energy absorption of a Y unit reads U ¼E h þE f ¼12M 0πC þð8βÀ4αÞM 0πH þ8γM 0H 2ð9ÞIn the quasi-static compression,the mean work of crushing force can be related to the total energy absorption by the Y unit.Here the effective crushing distance is not considered.So the power displacement is 2H :P m 2H ¼Uð10ÞSubstituting Eq.(9)into Eq.(10),P m can be determined by the following equation:P m ¼M 06πC H þ4γHt þð4βÀ2αÞπ ð11ÞWierzbicki pointed out that the independent parameters of thecollapse mode would minimize the crushing force [9].In other words,the differential equation dP m =dH is equal to 0,i.e.dP m ¼M 0À6πC H þ4γ¼0ð12ÞSoH ¼ffiffiffiffiffiffiffiffiffiffiffiffiffiffi32γπCt s ð13ÞP m can be obtained by substituting Eq.(13)into Eq.(11).And the mean crushing strength s m is equal to P m divided by the area A of Y unit section,that iss m ¼P m A ¼13ffiffiffiffiffiffiffiffiffi2γπp s 0t C 3=2þð2βÀαÞ6ffiffiffi3p s 0t C 2ð14Þwhere C ¼D =2,so s m ¼43ffiffiffiffiffiffiγπp s 0t D 3=2þ2ð2βÀαÞ3ffiffiffi3p s 0tD 2ð15ÞIt can be seen from Eq.(15)that the mean crushing strength has apower law relationship with t /D .De fining C 1¼4=3ffiffiffiffiffiffiγπp ,which is determined by adhesive failure length and plastic flowing area,and C 2¼2ð2βÀαÞ=ð3ffiffiffi3p Þwhich is determined by adhesive failure length and triangular pane,Eq.(15)can be simpli fied as follows:s m ¼C 1s 0t 32þC 2s 0t2ð16ÞEq.(16)shows that the mean crushing strength is dependent onmaterial plastic flowing stress;foil thickness and cell size of honeycomb.The first item of Eq.(16)is related to the plastic flowing area and plastic hinge line which is the function of cell edge length D ,while the second item related to the plastic hinge line which is determined by half-wavelength.3.3.Parameters calibrationThe coef ficients associated with the equations to calculate the mean crushing strength are generally obtained directly by assuming that they are constants related to half-wavelength H [8–13,15,16,20].However,such simpli fication is not justi fied in the current case because the effect of adhesive failure for deformation behavior much be taken into account.Here,the curve fitting technique is adopted to determine the parameters (C 1and C 2)in Eq.(16),based on the FE simulation results.3.4.FE model of a Y-section in the aluminum honeycombFE models have been extensively used to investigate the properties of honeycombs [4,19,21–24].A Y-section structure (Fig.8)was modeled using shell elements with the average mesh size of 0.3mm to obtain the mean crushing strength under quasi-static compression loading condition.The Y-section structure is placed between two rigid plates:one plate is fixed while the other is movable at a constant velocity in the axial direction,as shown in Fig.8.3.5.Model validationThe experimental data shown earlier is used to verify the FE models.The deformation patterns obtained from the model predictions and tests are compared in Fig.9.A reasonable agree-ment can be observed.Fig.10shows the computationally pre-dicted and experimentally measured nominal stress –strain responses.The results indicate a reasonable agreement.However,the fluctuation of the simulation nominal stress-strain curves is larger than that of the experimental data probably due to the excessive symmetric boundary constraints used for the Y-cell [19].The symmetric boundary conditions on three non-intersecting edges make the single cell-model actually equivalent to a honey-comb specimen consisting of repeated cells.Each fluctuation represents one fold formation of the cell wall in the simulation.However,in reality,the neighboring cells of a large sample could be interacted while forming the folds to reach their local peak value at different instants in the experiments,which results in smoother macroscopic experimental curves [19,25].The mean crushing strength calculated from the computational and experi-mental results are shown in Table 2.It is revealed that the predicted mean crushing strengths of Type 1,Type 2and Type 3samples are 2.1%higher,2.3%lower and 1.3%higher than thoseofFig.8.The FE model of a Y-section structure.Z.Bai et al./Thin-Walled Structures 80(2014)38–4542test data,respectively.Given the small discrepancy,the Y-section structure model can be deemed as validated.3.6.Coef ficients determinationUsing the validated model,a series of simulations were performed with different combinations of base materials and cell size to determine the unknown parameters in the newly devel-oped theoretical model.The model con figurations are listed in Tables 3and 4.In this study,a total number of 98cases were simulated.The loading and boundary conditions were consistent with those described in the above section.The effective plastic flowing stresses of Al5052and Al3003are 231.8MPa and 146.1MPa respectively [26–28].MATLAB (Version 7.0,The Mathworks,Inc.Natick,MA,U.S.A.)was employed forcurve fitting based on the data shown in Tables 3and 4to determine the values of coef ficients C 1and C 2using Eq.(16).As the result,C 1¼2:686,C 2¼7:946ðR 2¼0:9923Þwhich rely on simulation results obtained from the FE analysis and depend on the type of base material plastic flowing stress and the dimensions of the basic element.Then the new theoretical model to predict the mean crushing stress of hexagonal honeycomb is obtained,and it is written in Eq.(17).s m ¼2:686s 0t D 3=2þ7:946s 0tD 2ð17Þparison with similar studiesIn this section,the new formula developed was compared withthe theories by McFarland,Wierzbicki and Zarei Mahmoudabadi.The semi-empirical formula proposed by Mc Farland [5]has the following form:s m ¼s yst S 24:75K þ28:628 þ1:155q ys tS ð18Þwhere K is the ratio of 1/4wavelength and wall length D ;q ys is theshear yield stress of the base material;S is the cell minor diameter which is the distance between the two parallel cell edges;and t is the cell wallthickness.Fig.9.Deformation patterns obtained from the model prediction and test.0.00.51.01.52.02.5N o m i n a l S t r e s s (M P a )Nominal Strainputationally predicted and experimentally measured nominal stress –strain curves.Table 2The mean crushing strength of Type 1and 2aluminum honeycombs obtained based on the simulations and tests.Type numberMean crushing strength (MPa)Computationally predictedExperimentally measured Type 10.2400.245Type 2 1.184 1.157Type 30.3980.403Z.Bai et al./Thin-Walled Structures 80(2014)38–4543Zarei Mahmoudabadi and Sadighi [14–16]improved Wierzbicki's super folding element theory which under-estimates the mean crushing strength by considering the curvature effect.Fig.11shows the cylindrical curvature effect on the folding mode of a cell wall [15].The model by Zarei Mahmoudabadi is given below in the following equation:s m ¼s 0h 2ffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi3ðH Àb =2ÞS 2q 16b HI 1ðφ0Þþ6D αþ4H 2I 3ðφ0Þ"#ð19Þwhere b and αare the cylindrical curvature radius and angle,respectively,which can be estimated by experiment data.I 1ðφ0Þand I 3ðφ0Þare the integral expressions taken in the literatures [9,15].Mechanical properties data of five types of hexagonal alumi-num honeycomb with the same base material (Al5052H39)by HEXCEL CO.(Table 5)were employed as the baseline to compare the performance of the theoretical predictions.The material parameters of base material Al5052are documented in [26–29].The theoretical predictions of the mean crushing strength of the honeycomb based on the aforementioned models and the data in Table 4are compared in Fig.12.The least squares error methodwas applied to analyze the pared to the experi-mental data,the errors induced by the predictions based on the aforementioned models and new theoretical solutions are 19.99%(McFarland),14.89%(Wierzbicki), 2.65%(Zarei Mahmoudabadi)and 3.05%,respectively.The result of comparison shows that the accuracy of the improved model by Zarei Mahmoudabadi is slightly better than that of new formula,but the performance of the new formula is clearly higher than that of McFarland semi-empirical formula and Wierzbicki formula.In the current model development,the main factors inducing the error are as follows:(1)The base material was simpli fied as the ideal rigid-plastic inthe modeling.(2)Folding pattern was assumed to have a regular shape.Inreality,due to the imperfection of the base material,random folding pattern may occur.(3)FE model may introduce some errors.5.ConclusionThis paper presents a combined experimental,theoretical and computational study on the mean crushing strength of hexagonal multi-cell thin-walled structures.In the experimental study,three types of aluminum honey-comb samples were tested under the quasi-static compression loading conditions using an Instron machine.The nominal stress –strain curves were then obtained.The deformation pattern of the aluminum honeycomb was examined and used to guide the theoretical analysis.In the theoretical modeling,a new formula was derived based on an existing analytical model to predict the mean crushing strength of multi-cell thin-walled structures.The new formula shows that the mean crushing strength is a function of plastic flowing stress of the base material;foil thickness and cell size of honeycomb:numerical simulations were used to determine the unknown parameters in the new theoretical model,A largeTable 3Simulation matrix and results s m ðMPa Þof the honeycombs with base materialAl3003and different cell sizes.t (mm)D (mm)4.044.70 6.507.509.0011.0312.500.0400.4670.4100.23600.1830.1940.0970.0810.0520.7690.58090.34190.2720.2090.1560.1270.0600.9580.8740.4310.3590.2600.1800.1530.070 1.2110.9550.5990.4370.3340.2450.2030.080 1.539 1.2230.6930.5450.4260.2920.2490.100 2.153 1.676 1.0380.8470.5790.4220.3420.1202.9522.3731.3271.0100.8160.5890.456Table 4Simulation matrix and results s m ðMPa Þof the honeycombs with base materialAl5052and different cell sizes.t (mm)D (mm)4.044.70 6.507.509.0011.0312.500.0400.8080.6430.3850.3200.2280.1720.1440.052 1.3550.9440.5740.4610.3330.2500.2000.060 1.525 1.2720.7910.5860.4460.3120.2710.070 2.121 1.6670.9050.7350.5830.3980.3280.080 2.567 1.812 1.2140.9990.6310.4850.4100.100 3.770 2.9650.988 1.2791 1.0260.7480.6100.1205.0573.8522.3231.8611.3250.9510.730Fig.11.Cylindrical curvature effect on the folding mode of a cell wall:(a)the first type of folding mode in [14];and (b)the second folding mode in [15].Table 5Mechanical properties data of five types of hexagonal aluminum honeycomb with the same base material (Al5052H39)by HEXCEL CO.[18].Type number Cell size,S –alloy –foilthickness,t (mm)Nominal density (kg/m 3)Mean crushing strength (MPa)HEXCEL-1 3.969–5052–0.017841.60.621HEXCEL-2 4.763–5052–0.038170.4 1.724HEXCEL-3 6.350–5052–0.063583.2 2.310HEXCEL-49.525–5052–0.025425.60.276HEXCEL-59.525–5052–0.0508480.828Z.Bai et al./Thin-Walled Structures 80(2014)38–4544number of Y-section structure FE models with different con figura-tions were built and veri fied with the experimental data.The simulation results were then fitted with the newly developed formula to obtain the values of the coef ficients needed.Finally the predictions based on the new model and the other existing models were compared.The results show that the performance of new formula is much better than that of the theories by McFarland and Wierzbicki,but slightly worse than that by Zarei Mahmoudabadi.However the expression form of the new model is much simpler than the Zarei Mahmoudabadi model and requires much less parameters.Therefore,it is more suited in the practical structural design.It should be noted that the present model can only be used to predict the strength of hexagonal multi-cell thin walled structures with certain base material (Al5052and Al3003)and relatively small range of cell dimensions (D ¼4.04–12.5mm;t ¼0.04–0.12mm).This limitation is due to the fact that the two parameters C 1and C 2were determined based on the simulation results of the honeycombs with the aforementioned materials and dimensions.If the data are available for wider range parameters,the model can be extended to predict the behavior of more types of structures.AcknowledgmentsThis study is supported by the National Natural Science Foundation of China,under the Grant nos.51105137and 11172099.References[1]Zarei H,Kröger M.Optimum honeycomb filled crash absorber design.MaterDes 2008;29:193–204.[2]Lin Q,Nie H,Ren J,Chen JB.Investigation on design and reliability analysis of anew deployable and lockable mechanism.Acta Astronaut 2012;73:183–92.[3]Zhu F,Lu G,Ruan D,Shu D.Tearing of metallic sandwich panels subjected toair shock loading.Struct Eng Mech 2009;32:351–70.[4]Zhu F,Zhao L,Lu G,Gad E.A numerical simulation of the blast impact ofsquare metallic sandwich panels.Int J Impact Eng 2009;36:687–99.[5]Mc Farland R.Hexagonal cell structures under post-buckling axial load.AIAA J1963;1:1380–5.[6]Gibson L,Ashby M,Schajer G,Robertson C.The mechanics of two-dimensionalcellular materials.Proc R Soc Lond A:Math Phys Sci 1982;382:25–42.[7]Gibson LJ,Ashby MF.Cellular solids:structure and properties.Cambridgeuniversity press;1999.[8]De Oliveira JG,Wierzbicki T.Crushing analysis of rotationally symmetricplastic shells.J Strain Anal Eng Des 1982;17:229–36.[9]Wierzbicki T.Crushing analysis of metal honeycombs.Int J Impact Eng1983;1:157–74.[10]Wierzbicki T,Abramowicz W.On the crushing mechanics of thin-walledstructures.J Appl Mech 1983;50:727–34.[11]Wu E,Jiang W-S.Axial crush of metallic honeycombs.Int J Impact Eng1997;19:439–56.[12]Chen W,Wierzbicki T.Relative merits of single-cell,multi-cell and foam-filledthin-walled structures in energy absorption.Thin-Walled Struct 2001;39:287–306.[13]Zhang X,Cheng G,Zhang H.Theoretical prediction and numerical simulationof multi-cell square thin-walled structures.Thin-Walled Struct 2006;44:1185–91.[14]Zarei Mahmoudabadi M,Sadighi MA.Study on metal hexagonal honeycombcrushing under quasi-static loading.World Acad Sci,Eng Technol 2009;53:671–5.[15]Zarei Mahmoudabadi M,Sadighi M.A theoretical and experimental study onmetal hexagonal honeycomb crushing under quasi-static and low velocity impact loading.Mater Sci Eng:A 2011;528:4958–66.[16]Zarei Mahmoudabadi M,Sadighi M.A study on the static and dynamic loadingof the foam filled metal hexagonal honeycomb –theoretical and experimental.Mater Sci Eng:A 2011;530:333–43.[17]Bitzer T.Honeycomb technology:materials,design,manufacturing,applica-tions and testing.London:Chapman &Hall;1997.[18]Hexcel Corporation.HexWebTM honeycomb attributes and properties;November 1999.〈/Resources/DataSheets/Brochure-Da ta-Sheets/Honeycomb_Attributes_and_Properties.pdf 〉.[19]Hou B,Zhao H,Pattofatto S,Liu JG,Li YL.Inertia effects on the progressivecrushing of aluminium honeycombs under impact loading.Int J Solids Struct 2012;49:2754–62.[20]Liaghat G,Alavinia A.A comment on the axial crush of metallic honeycombsby Wu and Jiang.Int J Impact Eng 2003;28:1143–6.[21]Aktay L,Johnson AF,Kroplin BH.Numerical modelling of honeycomb corecrush behaviour.Eng Fract Mech 2008;75:2616–30.[22]Ali M,Qamhiyah A,Flugrad D,Shakoor M.Theoretical and finite elementstudy of a compact energy absorber.Adv Eng Softw 2008;39:95–106.[23]Shen CJ,Lu G,Yu TX.Dynamic behavior of graded honeycombs –a finiteelement pos Struct 2013;98:282–93.[24]Sun DQ,Zhang WH,Wei YB.Mean out-of-plane dynamic plateau stresses ofhexagonal honeycomb cores under impact pos Struct 2010;92:2609–21.[25]Wilbert A,Jang W-Y,Kyriakides S,Floccari J.Buckling and progressive crushingof laterally loaded honeycomb.Int J Solids Struct 2011;48:803–16.[26]Boyer HE,Gall TL,editors.Metals handbook.OH:Materials Park;1985.[27]Metals Handbook.Properties and selection:nonferrous alloys and special-purpose materials,vol.2,10th ed.:ASM International;1990.[28]Holt JMT,editor.Structural alloys handbook.West Lafayette,IN:CINDAS/Purdue University;1996.[29]Chawla A,Mukherjee S,Kumar D,Nakatani T,Ueno M.Prediction of crushingbehaviour of honeycomb structures.Int J Crashworth 2003;8:229–35.0.511.522.5312345M e a n c r u s h i n g s t r e s s (M P a )Typer of aluminum honeycombMcFarland [5]Wierzbicki [9]Zarei Mahmoudabadi [14-16]New ModelHexcel Experiment [18]parison of the theoretical model predictions and experimental data [18]on the mean crushing strength of aluminum honeycomb.Z.Bai et al./Thin-Walled Structures 80(2014)38–4545。

ANISOTROPIES IN THE COSMIC MICROW A VEBACKGROUND:THEORYSCOTT DODELSONNASA/Fermilab Astrophysics Center,P.O.Box500,Batavia,IL60510,USA Anisotropies in the Cosmic Microwave Background(CMB)contain a wealth of in-formation about the past history of the universe and the present values of cosmo-logical parameters.I ouline some of the theoretical advances of the last few years.In particular,I emphasize that for a wide class of cosmological models,theoristscan accurately calculate the spectrum to better than a percent.The specturm ofanisotropies today is directly related to the pattern of inhomogeneities present atthe time of recombination.This recognition leads to a powerful argument thatwill enable us to distinguish inflationary models from other models of structureformation.If the inflationary models turn out to be correct,the free parametersin these models will be determined to unprecedented accuracy by the upcomingsatellite missions.1HistoryThe Texas Symposium on Relativistic Astrophysics was held in Chicago ten years ago in1986.David Wilkinson spoke about the cosmic microwave back-ground.He undoubtedly made the point that the CMB provides us with some of the best evidence for the Big Bang.There was no evidence(and there still is no evidence)for any deviations from a black-body spectrum.And this is one of the primary predictions of the Big Bang.Wilkinson devoted most of his talk to searches for anisotropies in the CMB. The fact that the CMB temperature is the same in all directions indicates that the universe was very smooth early in its history.However,cosmologists generally work within the framework of gravitational instability which says that small inhomogeneities early on grew via gravity into the large structures we see today.Thus,the CMB should not be perfectly isotropic;it should carry some imprint of those small,early inhomogeneities.Wilkinson compiled the upper limits on anisotropies from the experiments of the time.This compilation is reproduced in Figure1,where I have taken the liberty of slightly changing his notation.In particular,it is convenient to expand the temperature on the sky in terms of spherical harmonicsT(θ,φ) T0=∞l=0lm=−la lm Y lm(θ,φ).(1)When we expand in this fashion,low l’s correspond to anisotropies on large an-1gular scales(the quadrupole is l=2)while large l’s correspond to anisotropies on small scales.The square of the coeffients of the Y lm’s are known as the C l’s.These are extremely useful things becuase they can be calculated by theorists and measured by observers.A given experiment at angular scale l measuresδT rms∼[l(l+1)C l/2π]1/2.The upper limits at the time correspond toδT rms∼50−200µK.Wilkinson was obviously aware of the fact that these upper limits were tantalizingly close to the levels of anisotropies predicted by many theories.He ended his talk by saying,“If the anisotropies are indeed just below current limits,as most of us feel they must be,the next few years should see thisfield turn from one of searching to one of studying.”2ExperimentsHow has thefield progressed since Wilkinson’s review?Figure1shows,along with Wilkinson’s compilation,a recent compilation1of all experiments in the last two years.Starting with the COBE2detection in1992at the largest angular scales,there have been dozens of detections on a wide range of angular scales.These detections,as Wilkinson’s quote makes clear,were anticipated based on typical models of structure formation.Although the details are not yet in,it is safe to say that gravitational instability theories predicted the level of anisotropies that are observed today.Another feature of the detections is just now becoming evident.As one moves from low l to high l(from large scales to small scales),one sees evidence of a gradual rise in the the amplitude of the anisotropies.We will see shortly that this too is a prediction of some of the more popular models of structure formation.How will the situation look ten years from now when results from the current crop of balloon-borne and ground-based experiments have come in, and the two satellite experiments(MAP3and PLANCK4)will have made all-sky maps?Figure1shows the expected error bars in the year2006.There are several ways to represent the knowledge we will have at the time.First, it is important to note that,today,experiments are sensitive to a range of l’s:thus,C500for example is not measured by a given small scale experiment. Rather,each experiment measures a signal integrated over a wide range of l. This range is depicted by the horizontal error bars in Figure1.In the future, we can continue to smooth over the l’s in this fashion.Then the errors will be as shown on the far right in Figure1.In order to see them on this graph, I have blown them up by a factor of100!We will also have the ability by then,though,to determine each individual C l.The expected errors on C5002Figure 1:Observations of the CMB spectrum.Upper limits are those compiled by Wilkinson in 1986.Detections were compiled by Lineweaver based on results within the last two years.Anticipated error bars from satellites are also shown.3are shown in Figure1.Either way you look at it,we will have an extraordinary amount of information in ten years.For the experimentalists,at least,it is clear that Wilkinson’s prediction has come true.Thefield really has moved from searching to studying.3TheoryCMB theorists have also been very active over the last few years.First of all,for a wide range of models,we are confident that we can calculate5the anisotropy spectra–the C l’s–to an accuracy of better than a percent.Figure2shows the results of seven different groups who independently calculated the C l’s for a given model.This graph was made about two years ago,and the agreement has only gotten better since then.Not only can we calculate accurately,but we can also calculate quickly.Thanks to Seljak and Zaldariagga6,in the time it has taken me to write this paragraph,we could have run offanother set of C l’s.We also have made great strides in the last few years understanding the bumps and wiggles in the theoretical curves.To understand the structure of these anisotropies,we need to review the thermal history of the universe.Recall that,early in the history of the universe,the temperature of the cosmic gas was very high.So,anytime a free electron and proton came together to form a hydrogen atom,a high energy photon immediately destroyed it.There was essentially no neutral hydrogen early on.This situation changed dramatically when the temperature dropped below1/3eV.After that time,there were not enough ionizing photons around.So almost all the free electrons and protons combined into neutral hydrogen.This had dramatic implications for the cosmic photons.As long as the electrons were free,they interacted with the photons via Compton scattering.After they combined into hydrogen,the photons travelled freely from the“surface of last scattering”to us today.So,for the purposes of the CMB,the universe is neatly divided into two epochs:Before Recombination when the photons and electrons behaved as a tightly coupled fluid and After Recombination when photons freestreamed.The mathematics of freestreaming is a little complicated,but the physics is completely trivial:it just requires us to trace the paths of free photons.So the physics behind the spectrum of anisotropies comes solely from the epoch Before Recombination.It pays to reiterate that Before Recombination,the photons and electrons acted as afluid.By this,I mean that it could be described by only its l=0 component(as opposed to all the multipole moments that are needed to de-scribe it today).This represents an immense simplification:instead of solving a infinite heirarchy of coupled differential equations for all the photon mo-40.0200.0400.0600.0800.01000.0Multipole Moments - l1.02.03.04.05.06.0l (l +1)C l / 2π (x 10-10)NOW !Figure 2:Seven different calculations of the CMB power spectrum for a model with adiabatic fluctuations.All lie within the shaded band which illustrates the minimum “cosmic variance”errors.ments,we need solve for only one of the moments.The forces acting on this moment,let’s call it δT ,are pressure and gravity.These forces act in oppo-site directions.Pressure tends to smooth out any inhomogeneities (i.e.drives δT to zero)while gravity produces inhomogeneities.It is not surprising then that acoustic oscillations are set up in the medium.In fact,Hu &Sugiyama 7have shown that this oscillation pattern is precisely the one imprinted in the C l spectrum of Figure 2.A quantitative analysis shows that there are two possible modes a that can be excited in this fluid.In particular,δT ( x ,η)= d 3ke i k · x A cos[kη/√+B sin[kη/√ (2)a Themodes look this simple only in the idealized case of zero baryons and pure matter domination.Accounting for baryons and other complications though does not alter the qualitative fact that there are two very distinct modes.5Figure3:The power for two different theories,one of which(adiabatic)excites the cosine mode of acoustic oscillations the other(isocurvature)the sine mode.Inflationary predictions typically look like the adiabatic spectrum here.Defect models should have some of the features of the isocurvature spectrum,but the calculations at present are not yet believable. whereηis conformal time.Again,not surprisingly,the C l spectrum today is radically different if the sine mode is excited than if the cosine mode is excited. Figure3shows that,as you would expect,the spectra are out of phase with each other.It is clear from the present data shown in Figure1that we will shortly be able to tell which of the two theoretical curves in Figure3is more accurate. That is,we will soon know whether the sine or the cosine mode were excited in the early universe.This is extremely important because we expect the two most popular mechanisms of structure formation–inflation and topological defects –to excite different modes.Let me walk through this argument which has recently been clearly elucidated by Hu&White8.Any theory which respects causality necessarily requires that there be no correlations on very large scales (scales that have not been in causal contact with each other).This is equivalent to a boundary condition onδT;namely that the Fourier transform vanishes at k=0.This means that only B in equation2can be non-zero.So topological6Figure4:The spectra for adiabatic models with different sets of parameters.Varied are the cosmological constant,neutrino mass,spectral index,and Hubble constant. defects,which of course obey causality,can be expected to excite the sine mode.Inflation is a theory which introduces correlations amongst scales that appear to be causally disconnected.Thus,inflationary models can,and most often do,excite the cosine mode.There are several caveats to the above argument.First of all,the pre-dictions for the adiabatic models depend on various cosmological parameters9: the slope of the primordial spectrum,the contribution from tensor modes, the Hubble constant,the baryon density,and several others.Thus the actual curves share some of the features of the curve labelled“Adiabatic”infigure3, but the predictions are by no means unique(seefigure4).Fortunately there are some robust features of these curves which hold up even after allowing many parameters to vary.The second caveat is that we simply do not know for sure that defect theories follow the general isocurvature model.There have been a few calculations of the spectrum in defect models10.As one who is actively at work on one such calculation,I think it is fair to say that we have not yet reached agreement.Assuming there are no major theoretical surprises,we can expect the ex-7periments over the next several years to pick out whether toppological defects or inflation are correct.Once that issue is settled,it remains to pin down the cosmological parameters which impact upon the spectrum.One might think that since there are so many free parameters,they cannot all be determined simultaneously.Recent work has shown that this is not true.Figure5shows an example11:we letfive parameters vary and show the error ellipses projected down onto a couple of two dimensional planes.The topfigure shows that it is quite possible that by the year2006,we will not be arguing about whether the Hubble constant is50or100,but rather whether it is50.5or50.0.The bottom figure shows that,in addition to the cosmological parameters,we should get a good handle on the inflationary parameters,thereby allowing us to distinguish amongst different inflationary models12.A number of groups13have varied even more parameters and all have reached the same general conclusion:the cosmological parameters will be pinned down to unprecedented accuracy by the satellite experiments.4ConclusionsAbout thirty years ago,Penzias and Wilson discovered the cosmic microwave radiation.This discovery convinced the vast majority of physicists that the Big Bang model was correct.In1992,the COBE satellite discovered anisotropies in the CMB.The existence and amplitude of these anisotropies were predicted by theories which relied on gravitational instability to form structure.It is perhaps too early to know for sure,but I would guess that COBE’s most enduring legacy will be its evidence that current models of structure formation are on the right track.A number of cosmologists are beginning to speculate about what we will have learned in ten years after the next generation of balloon and ground based experiments and after the MAP and PLANCK satellites haveflown.There is a very good chance these measurements will clearly distinguish between the two most popular models of structure formation:inflation and topological defects. Indeed,this could happen very soon.If a peak does indeed develop in the C l spectrum at around l∼200,this will be strong evidence for inflation.If the general picture of inflation is verified in this manner,the fun will begin. It will then be possible to determine many of the cosmological parameters to unprecedented accuracy.Further,the experiments will contain so much information that it will be possible to distinguish amongst different inflationary models.This opens a window to study physics at energies that are twelve orders of magnitude higher than those probed by the largest accelerators.Many people are fond of pointing out that something completely unex-8Figure 5:The estimated 95%contours for the MAP (larger ellipses in each case)and PLANCK (formerly COBRAS/SAMBA)satellites.In each case five variables –the am-plitude of the scalar and tensor perturbations,the spectral index of the scalars,the baryon density,and the Hubble constant –are allowed to vary.The ellipses are the projections of the five dimensional ellipses onto the (Hubble constant,Baryon density)plane and the (spectral index,tensor/scalar ratio)plane.The points and lines in the n −r plane correspond to thepredictions of different inflationary models.9pected and confusing may turn up,thereby upsetting the possibility of any such determinations.Of course this is possible.Unexpected and confusing discoveries have rocked cosmology for decades.The“man-bites-dog”story in cosmology though is the one in which the confusion ends;this may well happen within the next ten years.AcknowledgmentsThis work was supported in part by DOE and NASA grant NAG5–2788at Fermilab..References1.C.H.Lineweaver et al,astro-ph/9610133.2.G.P.Smoot et al,Astrophys.J.396,L1(1992).3.The MAP home page is /.4.The PLANCK(formerly COBRAS/SAMBA)home page ishttp://astro.estec.esa.nl/SA-general/Projects/Cobras/cobras.html.5.P.J.E.Peebles and J.T.Yu,Astrophys.J.162,815(1970);M.L.Wilsonand J.Silk,Astrophys.J.243,14(1981);J.R.Bond and G.Efstathiou,Astrophys.J.285,L45(1984).6.U.Seljak and M.Zaldariagga,Astrophys.J.469,437(1996).7.W.Hu and N.Sugiyama,Astrophys.J.444,489(1995); F.Atrio-Barandela and A.G.Doroshkevich,Astrophys.J.420,26(1994);P.Nasel’skij and I.Novikov,Astrophys.J.413,14(1993);U.Seljak,As-trophys.J.435,L87(1994);A.G.Doroshkevich.Ya.B.Zeldovich,andR.A.Sunyaev,Sov.Astron.22,523(1978);A.G.Doroshkevich,Sov.Astron.Lett.14,125(1988).8.W.Hu and M.White,Phys.Rev.Lett.77,1687(1996);N.G.Turok,astro-ph/9607109(1996).9.J.R.Bond,R.Crittenden,R.L.Davis,G.Efstathiou and P.J.Stein-hardt,Phys.Rev.Lett.72,13(1994).10.R.G.Crittenden and N.G.Turok,Phys.Rev.Lett.75,2642(1995);A.Albrecht,D.Coulson,P.Ferreira,and J.Magueijo,Phys.Rev.Lett.76,1413(1996).11.S.Dodelson,W.Kinney,and E.W.Kolb,(1997).12.L.Knox and M.S.Turner,Phys.Rev.Lett.73,3347(1994).13.L.Knox,Phys.Rev.D52,4307(1995);G.Jungman,M.Kamionkowski,A.Kosowsky,and D.N.Spergel,Phys.Rev.D54,1332(1996);S.Dodelson,E.I.Gates,and A.S.Stebbins,Astrophys.J.467,10(1996).10。

a r X i v :0805.2757v 1 [q u a n t -p h ] 18 M a y 2008Robustness of Adiabatic Quantum Computing Seth Lloyd,MIT Abstract:Adiabatic quantum computation for performing quantum computations such as Shor’s algorithm is protected against thermal errors by an energy gap of size O (1/n ),where n is the length of the computation to be performed.Adiabatic quantum computing is a novel form of quantum information processing that allows one to find solutions to NP-hard problems with at least a square root speed-up [1].In some cases,adiabatic quantum computing may afford an exponential speed-up over classical computation.It is known that adiabatic quantum computing is no stronger than conventional quantum computing,since a quantum computer can be used to simulate an adiabatic quantum computer.Aharonov et al.showed that adiabatic quantum computing is no weaker than conventional quantum computation [2].This paper presents novel modelsfor adiabatic quantum computation and shows that adiabatic quantum computation is intrinsically protected against thermal noise from the environment.Indeed,thermal noise can actually be used to ‘help’an adiabatic quantum computation along.A simple way to do adiabatic versions of ‘conventional’quantum computing is to use the Feynman pointer consisting of a line of qubits [3].The Feynman Hamiltonian isH =−n −1ℓ=0U ℓ⊗|ℓ+1 ℓ|+U †ℓ⊗|ℓ ℓ+1|,(1)where U ℓis the unitary operator for the ℓ’th gate and |ℓ is a state of the pointer where the ℓ’th qubit is 1and all the remaining qubits are 0.Clearly,H is local and each of itsterms acts on four qubits at once for two-qubit gates.If we consider the pointer to be a ‘unary’variable,then the each of the terms of H acts on two qubits and the unary pointer variable.Assume that the computation has been set up so that all qubits start in the state 0.The computation takes place and the answer is placed in an answer register.Now a long set of steps,say n/2takes place in which nothing happens.Then the computation is undone,returning the system to its initial state at the n−1’th step.The computational circuit then wraps around and begins again.The eigenstates of H then have the form|b,k =(1/√Now gradually turn on an H term while turning offthe H1term:the total Hamiltonian isηH0+(1−λ)H1+λH.Asλis turned on,the system goes to its new ground state |b=0,k=0 .It can be verified numerically and analytically[7]that the minimum energygap in this system occurs atλ=1:consequently,the minimum gap goes as1/n2.In fact, the energy gap due to the interaction between the H1and the H terms is just the energy gap of the simpler system consisting just of the chain qubits on their own,confined to the subspace in which exactly one qubit is1:that is,it is the energy gap of a qubit chain with Hamiltonian(1−λ)H1+λH′,where H′=− ℓ|ℓ+1 ℓ|+|ℓ ℓ+1|.This gap goes as1/n2.Accordingly,the amount of time required to perform the adiabatic passage is polynomial in n.When the adiabatic passage is complete,the energy gap of the H term on its own goes as1/n2from the cosine dependence of the eigenvalues of H:it is also just the energy gap of the simplified system in the previous paragraph.This implies that the adiabatic passage can accurately be performed in a time polynomial in n.Measuring the answer register now gives the answer to the computation with probability1/2.This is an alternative(and considerably simpler)demonstration to that of[2]that‘conventional’quantum computation can be performed efficiently in an adiabatic setting.An interesting feature of this procedure is that the adiabatic passage can be faulty and still work justfine:all of the energy eigenstates in the|b=0,k sector give the correct answer to the computation with probability1/2,for any k.The real issue is making sure we do not transition to the|b=0,k sector.But the Hamiltonians H1and H do not couple to this sector:so in fact,we can perform the passage non-adiabatically and still get the answer to the computation.For example,if we turn offthe H1Hamiltonian very rapidly and turn on the H Hamiltonian at the same time,the system is now in an equal superposition of all the|b=0,k eigenstates.If we wait for a time∝n2(corresponding to the inverse of the minimum separation between the eigenvalues of H),then the state of the system will be spread throughout the|b=0,k sector,and we can read out the answer with probability1/2.This method effectively relies on dispersion of the wavepacket tofind the answer.Since coherence of the pointer doesn’t matter,we can also apply a Hamiltonian to the pointer that tilts the energy landscape so that higher pointer values have lower energy.(E.g.,H pointer=− ℓℓE|ℓ ℓ|.)Starting the pointer offin the initial state above and letting it interact with a thermal environment will obtain the answer to the computationin time of O(n).Similarly,in the absence of an environment,starting the pointer offin a wavepacket with zero momentum at time0and letting it accelerate will get the answer to the computation in even shorter time.Clearly,this method is quite a robust way of performing quantum computation.Let us look more closely at the sources of this robustness.Ifηis big,then there is a separation of energy of O(η/n)between the|b=0,k sector—states which give the correct answer to the computation—and the|b=0,k sector—states which give the incorrect answer to the computation.This is because b=0,k|ηH0|b=0,k =η/n.This energy gap goesdown only linearly in the length of the computation and can be made much larger than the gap between the ground andfirst excited state by increasingη>>1.This second energy gap is very useful:it means that thermal excitations with an energy below the gap will not interfere with obtaining the proper answer.That is,this method is intrinsically error-tolerant in the face of thermal noise.Indeed,it is this O(η/n) gap that determines how rapidly the computation can take place rather than the O(1/n2) gap between the ground and excited states.Of course,the actual errors in a system that realizes the above scheme are likely to arise from variability in the applied Hamiltonians.The energy gap arguments for robustness only apply to the translational dynamics of the system(this is what makes the analysis of the system tractable in thefirst place).That is,errors that affect each Uℓon its own are not protected against:but these are the errors that cause the computation to come out wrong.Of course,one can always program the circuit to perform‘conventional’robust quantum computation to protect against such errors.One must be careful,however,that errors that entangle the pointer with the logical qubits do not contaminate other qubits: conventional robust quantum computation protocols will have to be modified to address this issue.Farhi et al.have recently exhibited error correcting codes for‘conventional’adiabatic quantum computation[1]that can protect against such computational errors[8].The use of error correcting codes to correct the variation in the Uℓmay well be overkill. In any system manufactured to implement adiabatic quantum computing,these errors in the Uℓare essentially deterministic:the Uℓcould in principle be measured and their varia-tion from their nominal values compensated for by tuning the local Hamiltonians.Because it involves no added redundancy,such an approach is potentially more efficient than the use of quantum error correcting codes.Exactly how to detect and correct variations in the Uℓwill depend on the techniques(e.g.,quantum dots or superconducting systems)used toconstruct adiabatic quantum circuits.It is also interesting to note that performing quantum computation adiabatically is intrinsically more energy efficient than performing a sequence of quantum logic gates via the application of a series of external pulses.The external pulses must be accurately clocked and shaped,which requires large amounts of energy.In the schemes investigated here,the internal dynamics of the computer insure that quantum logic operations are performed in the proper order,so no clocking or external pulses need be applied.The adiabatic technique also avoids the Anderson localization problem raised by Landauer.The above construction requires an external pointer and four qubit interactions.One can also set up a pointerless model that requires only pairwise interactions between spin1/2 particles(compare the following method with the method proposed in reference[9]).Let each qubit in the computational circuit correspond to a particle with two internal states. Let each wire in the circuit correspond to a mode that can be occupied by a particle.The ℓ’th quantum logic gate then corresponds to an operator˜Hℓ=Aℓ+A†,where Aℓis anℓoperator that takes two particles from the two input modes and moves them to the output modes while performing a quantum logic operation on their two qubits.That is,Aℓ=a†out1a in1a†out2a in2⊗Uℓ.(5) Note that Aℓacts only when both input modes are occupied by a qubit-carrying particle. If we use the Hamiltonian˜H= ℓ˜Hℓin place of H in the construction above,the ground state of this Hamiltonian is a superposition of states in which the computation is at various stages of completion.Just as above,measurement on the ground state will reveal the answer to the computation with probability1/2.Note that even though the Hamiltonian in equation(5)involves a product of opera-tors on four degrees of freedom(the internal degrees of freedom of the particles together with their positions),it is nonetheless a physically reasonable local Hamiltonian involv-ing pairwise interactions between spin-1/2particles.To simulate its operation using an array of qubits as in[2]would require four qubit interactions,as in the pointer model dis-cussed above.This point is raised here because of the emphasis in the quantum computing literature on reducing computation to pairwise interactions between qubits.Pairwise in-teractions between particles orfields–i.e.,the sort of interactions found in nature–may correspond to interactions between more than two qubits.Without further restrictions on the form of the quantum logic circuit,evaluating the energy gap in this particle model is difficult,even for thefinal Hamiltonian1−˜H.But wecan always set up the computational circuit in a way that allows the adiabatic passage to be mapped directly on the Feynman pointer model above.The method is straightforward: order the quantum logic gates as above.Now insert additional quantum logic gates between each consecutive pair of gates in the original circuit.The additional gate inserted between theℓ’th andℓ+1’th quantum logic gates couples one output qubit of theℓ’th quantum logic gate with one input qubit of theℓ+1’th gate,and performs a trivial operation U=1on the internal qubits of these gates.The purpose of these gates is ensure that the quantum logic operations are performed in the proper sequence.Effectively,one of the qubits from theℓ’th gate must‘tag’one of the qubits from theℓ+1’th gate before theℓ+1’th gate can be implemented.Accordingly,we call this trick,a‘tag-team’quantum circuit.Tag-team quantum circuits are unitarily equivalent to the Feynman pointer model with an extra,identity gate inserted between each of the original quantum logic gates. Accordingly,the spectral gap for tag-team quantum circuits goes as1/n2and the quantum computation can be performed in time O(poly(n)).Just as for the pointer version of adiabatic quantum computing,the important spectral gap for tag-team adiabatic quantum computation is not the minimum gap,but rather the gap of size O(η/n)between the ground-state manifold of‘correct’states and the next higher manifold of‘incorrect’states. Once again,the existence of this gap is a powerful protection against errors in adiabatic quantum computation.The methods described above represent an alternative derivation of the fact that adiabatic quantum information processing can efficiently perform conventional quantum computation.The relative simplicity of the derivation from the original Feymnan Hamil-tonian[3]allows an analysis of the robustness of the scheme against thermal excitations. Adiabatic implementations of quantum computation are robust against thermal noise at temperatures below the appropriate energy gap.The appropriate energy gap is not the minimum gap,which scales as n2,but the gap between the lowest sector of eigenstates, which give the correct answer,and the next sector.This gap scales asη/n,whereηis an energy parameter that is within the control of the experimentalist.Acknowledgements:This work was supported by ARDA,ARO,DARPA,CMI,and the W.M.Keck foundation.The author would like to thank D.Gottesman and D.Nagaj for helpful conversations.References:[1]E.Farhi,J.Goldstone,S.Gutmann,Science292,472(2001).[2]D.Aharonov,W.van Dam,J.Kempe,ndau,S.Lloyd,O.Regev,‘Adiabatic Quantum Computation is Equivalent to Standard Quantum Computation,’Proceedings of the45th Annual IEEE Symposium on Foundations of Computer Science(FOCS’04), 42-51(2004);quant-ph/0405098.[3]R.Feynman,Found Phys.16,507-531(1986).[4]E.Farhi,S.Gutmann,Phys.Rev.A58,915-928(1998).[5]D.Aharonov,A.Ambainis,J.Kempe,U.Vazirani,Proceeding of the thirty-third annual ACM symposium on Theory of Computing,pp.50-59(2001).[6]ndauer,Phil.Trans.:Phys.Sci.and Eng.,353,367-376(1995).[7]P.Deift,M.B.Ruskai,W.Spitzer,arXiv:quant-ph/0605156.[8]S.P.Jordan,E.Farhi,P.Shor,Phys Rev A74,052322(2006);quant-ph/0512170.[9]A.Mizel,D.A.Lidar,M.Mitchell,Phys.Rev.Lett.99,070502(2007)arXiv:quant-ph/0609067.。

原子与分子物理学报JOURNAL OF ATOMIC AND MOLECULAR PHYSICS第38卷第1期2021年2月Vol. 38 No. 1Feb. 2021里德伯钠原子经绝热快速通道的布居数跃迁相干控制蒋利娟，王金龙，赫丙玲，田时振(新乡学院物理与电子工程学院，新乡453003)摘要：运用含时多态展开方法和B-样条函数研究了 O 啾频率微波场中里德伯钠原子的量子态之间的布居数迁移.计算了里德伯钠原子n = 70 -77的开普勒频率.计算了在不同的微波场中六个态的布居数从n = 70到n = 77随时间的迁移，布居数跃迁到最终态n = 77达到了 98%，这可以通过连续的单光子跃迁来实现.结果表明，通过优化微波脉冲参数可以实现从低态到较高态的布居跃迁的相干控制.关键词：B-样条；含时多态展开方法；相干控制；微波脉冲中图分类号：O562. 1 文献标识码：A DOI ： 10.19855/j.l000-0364.2021.011004Numerical exploration of population transfer ofRydberg sodium atom via adiabatic passageJIANG Li-Juan , WANG Jin-Long , HE Bing-Ling , TIAN Shi-Zhen(Department of Physics and Electronic Engineering , Xinxiang University , Xinxiang 453003 , China)Abstract : We demonstrate theoretically that Rydberg sodium atoms can be transferred to states of higher princi ­pal quantum number by exposing them to specially designed frequency 一 chirped Gaussian microwave pulses. Thepopulation has been transferred from n = 70 to n = 77 states of sodium atoms with efficiency of more than 98% , which is achieved by means of the sequential single 一 photon A n 二一 1 transitions. The results show that the co ­herent control of the population transfer from the lower n to the higher n states can be accomplished by optimiza ­tion of the microwave field parameters with the time 一 dependent multilevel approach (TDMA ).Key words : B - spline ; Time - dependent multilevel approach (TDMA ) ； Coherent control ; Microwave pulse1引言近年来，激光脉冲的设计在量子态之间有效 实现且产生较大的布居迁移已成为许多理论和实验研究的主题［1-11］.该课题涉及许多应用，包括 光谱学、碰撞动力学和化学反应的控制.可调谐 染料激光器允许高里德伯原子的n 或2态达到有效的布居• Meerson 和Friedlang ［12］发现，运用微 波脉冲，最初在开普勒频率并碉啾到较低频率,将原子转移到更高的n 态，这导致在较低的微波场电离.Bensky 等［13〕和Wesdorp 等［14〕建议在另 一个方向碉啾频率可能更有意思，他们证明了运用半周期脉冲将电子重组引入高位里德伯状态是可以实现的，该技术可能是产生反氢一种方法［15］-最近，通过频率或振幅碉啾脉冲诱导里德伯原子的相干布居转移被运用到物理和化学过程中量子态之间的控制［16'17］ •结果表明，激光脉冲参数对布居跃迁几率起决定性作用，因此研究激光脉冲参数对布居跃迁几率的提高具有重要意义.绝热快速通道(Adiabatic rapid passage )可以实现布居数在多态系统中的100%跃迁，是布居数从一个量子态到另一个量子态迁移的最有效方法⑴，9L 运用碉啾辐射场可以使布居迁移，通过多个绝热快速通道可以使布居数达到有效的迁移.收稿日期：2019-12-14基金项目：国家自然科学基金(U1704132)；河南省高校科技创新人才项目(14HASTIT044)；河南省教育厅科学技术研究(20B140013)；河南省高等学校重点科研项目(21A430031)；新乡学院“理论物理”重点学科建设(2018119)作者简介：蒋利娟(1972—)，女，河南滑县人，教授，博士，研究方向为原子与分子物理.E-mail ： *******************第38卷原子与分子物理学报第1期最近，运用单脉冲的绝热快速通道把布居数囚禁在多级梯形原子中已被广泛讨论•Lambert等人旦]证明了通过运用频率碉啾微波脉冲，可以把里德伯原子转移到较低的量子状态•在他们的实验中，运用碉啾7.8-11.8GHz之间的脉冲，在不同态之间通过双光子共振，里德伯锂原子的布居数从n=75迁移到n=66.Djotyan等人卩1]设计了一种方法，运用单频率碉啾激光脉冲，没有考虑原子的激发，使八原子的两个基态之间进行布居数跃迁•当激光脉冲包络频率光谱(没有碉啾)的宽度较小时，脉冲的Rabi频率峰值大于子的两基的频率，可实数跃迁.在他们的分析研究中，两个过渡偶极子元素假设是相等的，即丨九I=1匾3I[21]-但对于一个真实的系统，这并不总是令人满意.在本文中，我们运用设定的微波脉冲计算了里德伯钠原子的从n=70到77的布居数跃迁.为了获得更精确的数值结果，我们运用了许多耦合状态来代替所选择涉及到的态.通过解定态薛定谭方程，我们期望获得多能级系统的完全布居迁移的条件.2理论我们运用含时多态展开方法，该方法由Zhang等人】”〕提出•运用高斯微波脉冲[23]为：E(t)=E0e-21n2((!-!0)/!»)2cos W(i)i(1)其中,3(t)=盹+/31,®0为激光脉冲的初始频率, 0是碉啾率，E o、0、t”分别为微波脉冲的振幅、脉冲周期中心、脉冲周期宽度•钠原子的最外层电子的哈密顿量为：H=H o+zE(t)(2)式中H o是忽略原子自旋-轨道相互作用的零场哈密，方(外，文用原子单位)为：[-2!2+卩()]0()=E<p(r)，(3)在这里，卩()是碱金属原子价电子感受到原子实的势.包括渗透和极化的影响，可以写成下面的式子4]：V=-[(Z-1)exp(-5厂)+勺厂exp(-巾厂)+1] r(4)对于钠原子，Z二11,q=7.902，勺=23.51，5=2.688.零场下钠原子的波函数认(为：处“(r)=匸山(0，)其中n，/,m为主量子数、角量子数和磁量子数. Y lm(0")为球谐函数，人(r)为径向波函数•径向波函数人(r)是径向约化薛定谭方程的解：[-养+1{1^)+卩()]人(r)=E"人(r)(5)以B-样条函数为基函数，径向波函数人(r)可展开为：nU ni(r)=X(r)(6)i=1运用U ni(r)(r=0)=0，在这里B，(r)[25-30]是%阶的第i个B样条•将(6)代入(5)式可得到一方：H i C=ESC(7)这里H是哈密顿矩阵，S是重叠矩阵，E和C是本征 值和.通方(7)我可获E 和U ni(r).钠原子在微波场中时的波函数可以写成：n&(r，t)=X叫(t)以-遍"⑻i=1在上面的等式中，«ii(t)表示迁移几率的幅度的扩展系数，e-iEiit包含共振发生时动量中的相位信息[22]，波函数&(r,t)可以用含时薛定谭方程:i°叮：，)=H0(r,t)(9)d t将方程(8)代入方程(9)，求解含时薛定谭方程,我们可以得到切(t)•那么外层电子在i态的概率可：i=12P,=X I5I(10)i=0运用上面的公式，我们可以获得态与态跃迁的几率，也可以直接观察到著名的Rabi振荡和其他一些有趣的特征.3结果和讨论首先，我们讨论了微波脉冲参数0,t0和t”对n=73p态的布居跃迁的影响.设计的微波脉冲是高斯微波脉冲：E(t)=E0e-21n2((t-t0)/t»)2co SW(t)t,其中®(t)= %-0t.我们假定最初始态为n= 73p，选择®0=21.0GHz，E0=2.0V/cm,t0= 1000ns，t”=1500ns，0.0004GHz/ns W0W 0.0055GHz/ns.所获得的原子状态之间的布居迁移如图1所示.我们可以看到布居跃迁几率与碉啾率有关•因此，为了控制布居跃迁，碉啾率必第38卷蒋利娟，等：里德伯钠原子经绝热快速通道的布居数跃迁相干控制第1期须被优化.然后,我们选择0=0.0015GHz/ns和脉冲中心参数400ns W t0W1600ns,其他参数保持不变.获得布居跃迁如图2所示.我们可以看到在其他参数不变的情况下，n=73p态的布居跃迁几率随着微波脉冲中心参数S的增加而增加，S 大约为1000ns时，布居跃迁几率达到最大值,然后逐渐降低，直到t0〉1600ns时，布居跃迁几率基本为零•因此，里德伯钠原子的跃迁几率与微波脉冲中心参数t0有关.在图3中，我们选取了0=0.0015GHz/ns,t0=1000ns以及其他参数和图1中的参数保持一致，研究脉冲宽度t w （800ns W t w W2600ns）对跃迁几率的影响.从图3可以看到，布居跃迁几率随着微波脉冲宽度的改变而改变，这也表明微波脉冲宽度对布居跃迁几率是非常重要的参数•因此，为了相‘控制布居跃迁，脉冲宽度和其他的微波参数一样必须被优化.uo匸e-ndodChirped rate/(GHz/ns)图1选择e0=21.0GHz,E0=2.0V/cm,t0= 1000ns,t w=1500ns,子的(n=73)的布居数随碉啾率0的变化图像•Fig.1Population of the intermediate state(n=73)of sodium atom as a function of the chirped rate0,the other parameters are e0=21.0GHz,E0=2.0V/cm,t0=1000ns,t w=1500ns.现在我们讨论里德伯钠原子从n=70到n= 77态随碉啾脉冲19.2GHz—14.4GHz的布居跃迁几率.在这个过程中，通过运用频率脉冲进行扫描，一系列连续的绝热通道形成，主量子数n 和n+1的A n=1之间的跃迁主要发生在开普勒频率m/2n二1/2n n3⑶】处.从n=70到n=77的开普勒频率列在表1中.在图4中，通过优化这些参数，我们选择®=21.0GHz,E°=2.0V/ cm,仔=0.0015GHz/ns,t0=1000ns,t w=1500 ns.在图4中，很明显的布居迁移可以观察到,86421aaao.0.uo匸e-ndod4006008001000120014001600t0/ns图2选择盹=21.0GHz,E0=2.0V/cm,0=0.0015GHz/ns,t w=1500ns,钠原子的中间态(n=73)的布居数随脉冲中心参数t。

a rXiv:h ep-ph/6614v212J un26THE SEARCH FOR DARK MATTER AXIONS P.SIKIVIE Theoretical Physics Division,CERN,CH-1211Gen`e ve 23,Switzerland and Department of Physics,University of Florida,Gainesville,FL 32611,USA Axions solve the Strong CP Problem and are a cold dark matter candidate.The combined constraints from accelerator searches,stellar evolution limits and cosmology suggest that the axion mass is in the range 3·10−3>m a >10−6eV.The lower bound can,however,be relaxed in a number of ways.I discuss the constraint on axion models from the absence of isocurvature perturbations.Dark matter axions can be searched for on Earth by stimulating their conversion to microwave photons in an electromagnetic cavity permeated by a magnetic fiing this technique,limits on the local halo density have been obtained by the Axion Dark Matter eXperiment.1IntroductionThe standard model of elementary particles has been an extraordinary breakthrough in our description of the physical world.It agrees with experiment and observation disconcertingly well and surely provides the foundation for all further progress in our field.It does however present us with a puzzle.Indeed the action density includes,in general,a termL stand mod =...+θg 2has no such U A(1)symmetry2.Second,thatθis cyclic,that is to say that physics atθis indistinguishable from physics atθ+2π.Third,that an overall phase in the quark mass matrix m q can be removed by a redefinition of the quarkfields only if,at the same time,θis shifted to θ−arg det m q.The combination of standard model parameters¯θ≡θ−arg det m q is independent of quarkfield redefinitions.Physics therefore depends onθsolely through¯θ.Since physics depends on¯θ,the value of¯θis determined by experiment.The term shown in Eq.(1)violates P and CP.This source of P and CP violation is incompatible with the experimental upper bound on the neutron electic dipole moment unless|¯θ|<10−10.A new improved upper limit on the neutron electric dipole moment(|d n|<3.0·10−26e cm)was reported at this conference by P.Geltenbort3.The puzzle aforementioned is why the value of¯θis so small.It is usually referred to as the“Strong CP Problem”.If there were only strong interactions,a zero value of¯θcould simply be a consequence of P and CP conservation.That would not be much of a puzzle.But there are also weak interactions and they,and therefore the standard model as a whole,violate P and CP.So these symmetries can not be invoked to set¯θ=0.More pointedly,P and CP violation are introduced in the standard model by letting the elements of the quark mass matrix m q be arbitrary complex numbers4.In that case,one expects arg det m q,and hence¯θ,to be a random angle.The puzzle is removed if the action density is insteadL stand mod+axion=...+132π2a(x)f a−det arg m q depends on the expectation value of a(x).That expectation value minimizes the effective potential.The Strong CP Problem is thensolved because the minimum of the QCD effective potential V(¯θ)occurs at¯θ=07.The weak interactions induce a small value for¯θ8,9,of order10−17,but this is consistent with experiment.The notion of Peccei-Quinn(PQ)symmetry may seem contrived.Why should there be a U(1)symmetry which is broken at the quantum level but which is exact at the classical level? However,the reason for PQ symmetry may be deeper than we know at present.String theory contains many examples of symmetries which are exact classically but which are broken by quantum anomalies,including PQ symmetry10,11,12.Withinfield theory,there are examples of theories with automatic PQ symmetry,i.e.where PQ symmetry is a consequence of just the particle content of the theory without adjustment of parameters to special values.Thefirst axion models had f a of order the weak interaction scale and it was thought that this was an unavoidable property of axion models.However,it was soon pointed out13,14that the value of f a is really arbitrary,that it is possible to construct axion models with any value of f a.A value of f a far from any previously known scale need not lead to a hierarchy problem because PQ symmetry can be broken by the condensates of a new technicolor-like interaction 15.The properties of the axion can be derived using the methods of current algebra16.Theaxion mass is given in terms of f a bym a≃6eV106GeVπa(x)f f =ig fm ffγ5f(5)where g f is a model-dependent coefficient of order one.In the KSVZ model the coupling toelectrons is zero at tree level.Models with this property are called’hadronic’.The axion has been searched for in many places but not found17.The resulting constraintsmay be summarized as follows.Axion masses larger than about50keV are ruled out by particlephysics experiments(beam dumps and rare decays)and nuclear physics experiments.The next range of axion masses,in decreasing order,is ruled out by stellar evolution arguments.Thelongevity of red giants rules out200keV>m a>0.5eV18,19in the case of hadronic axions,and200keV>m a>10−2eV20in the case of axions with a large coupling to electrons[g e=0(1)in Eq.5].The duration of the neutrino pulse from Supernova1987a rules out2 eV>m a>3·10−3eV21.Finally,there is a lower limit,m a>∼10−6eV,from cosmology which will be discussed in detail in the next section.This leaves open an“axion window”:3·10−3>m a>∼10−6eV.We will see,however,that the lower edge of this window(10−6eV) is much softer than its upper edge.2Axion cosmologyThe implications of the existence of an axion for the history of the early universe may be briefly described as follows.At a temperature of order f a,a phase transition occurs in which the U P Q(1)symmetry becomes spontaneously broken.This is called the PQ phase transition.At these temperatures,the non-perturbative QCD effects which produce the effective potential V(θ)turns on and the axion acquires mass.There is a critical time,defined by m a(t1)t1=1,when the axionfield starts to oscillate in response to the turn-on of the axion mass.The corresponding temperature T1≃1 GeV23.The initial amplitude of this oscillation is how far from zero the axionfield lies when the axion mass turns on.The axionfield oscillations do not dissipate into other forms of energy and hence contribute to the cosmological energy density today23.This contribution is called of‘vacuum realignment’.It is further described below.Note that the vacuum realignment contribution may be accidentally suppressed in case1if the homogenized axionfield happens to lie close to zero.In case2the axion strings radiate axions24,25from the time of the PQ transition till t1when the axion mass turns on.At t1each string becomes the boundary of N domain walls.If N=1,the network of walls bounded by strings is unstable26,27and decays away.If N>1there is a domain wall problem28because axion domain walls end up dominating the energy density, resulting in a universe very different from the one observed today.There is a way to avoid this problem by introducing an interaction which slightly lowers one of the N vacua with respect to the others.In that case,the lowest vacuum takes over after some time and the domain walls disappear.There is little room in parameter space for that to happen and we will not consider this possibility29further here.Henceforth,we assume N=1.In case2there are three contributions to the axion cosmological energy density.One contri-bution24,25,30,31,32,33,34is from axions that were radiated by axion strings before t1.A second contribution is from axions that were produced in the decay of walls bounded by strings after t131,35,36,29.A third contribution is from vacuum realignment23.Let me briefly indicate how the vacuum alignment contribution is evaluated.Before time t1, the axionfield did not oscillate even once.Soon after t1,the axion mass is assumed to change sufficiently slowly that the total number of axions in the oscillations of the axionfield is an adiabatic invariant.The average number density of axions at time t1isn a(t1)≃1t1(6)In Eq.(6),we used the fact that the axionfield a(x)is approximately homogeneous on the horizon scale t1.Wiggles in a(x)which entered the horizon long before t1have been red-shifted away37.We also used the fact that the initial departure of a(x)from the nearest minimum is of order f a.The axions of Eq.(6)are decoupled and non-relativistic.Assuming that the ratio of the axion number density to the entropy density is constant from time t1till today,onefinds 23,29Ωa≃11012GeV 7h 2(7) for the ratio of the axion energy density to the critical density for closing the universe.h is the present Hubble rate in units of100km/s.Mpc.The requirement that axions do not overclose the universe implies the constraint m a>∼6·10−6eV.The contribution from axion string decay has been debated over the years.The main issue is the energy spectrum of axions radiated by axion strings.Battye and Shellard32have carried out computer simulations of bent strings(i.e.of wiggles on otherwise straight strings)and have concluded that the contribution from string decay is approximately ten times larger than that from vacuum realignment,implying a bound on the axion mass approximately ten times more severe,say m a>∼6·10−5eV instead of m a>∼6·10−6eV.My collaborators and I have done simulations of bent strings31,of circular string loops31,34and non-circular string loops34.We conclude that the string decay contribution is of the same order of magnitude than that from vacuum realignment.Yamaguchi,Kawasaki and Yokoyama33have done computer simulations of a network of strings in an expanding universe,and concluded that the contribution from string decay is approximately three times that of vacuum realignment.The contribution from wall decay has been discussed in detail in ref.29.It is probably subdominant compared to the vacuum realignment and string decay constributions.It should be emphasized that there are many sources of uncertainty in the cosmological axion energy density aside from the uncertainty about the contribution from string decay.The axion energy density may be diluted by the entropy release from heavy particles which decouple before the QCD epoch but decay afterwards38,or by the entropy release associated with a first order QCD phase transition.On the other hand,if the QCD phase transition isfirst order 39,an abrupt change of the axion mass at the transition may increaseΩa.If inflation occurs with reheat temperature less than T P Q,there may be an accidental suppression ofΩa because the homogenized axionfield happens to lie close to a CP conserving minimum.Because theRHS of Eq.(7)is multiplied in this case by a factor of order the square of the initial vacuum misalignment angle a(t1)x.Recently,Kaplan and Zurek proposed a model 40in which the axion decay constant fa is time-dependent,the value f a(t1)during the QCD phase-transition being much smaller than the value f a today.This yields a suppression of the axion cosmological energy density by a factor(f a(t1)<δa( x,t)δa( x′,t)>e−i k·( x− x′)= H I k3,(8)(2π)3where x are comoving coordinates.Eq.(8)is often written in the shorthand notationδa=H Iα( x,t1)2(9)2t1whereα( x,t1)=a( x,t1)/f a is the local misalignment angle.Thefluctuations in the axionfieldproduce perturbations in the axion dark matter densityδn iso aa1=H I2(g00−1).Note that in case1the density perturbations in the cold axionfluid have bothadiabatic and isocurvature components.The adiabatic perturbations(δρad a4ρr=δTρCDM =δρiso aρCDM=H IΩCDM<0.31δρmρm is the amplitude of the primordial spectrum of matter perturbations.It is related to the amplitude of large scale(low multipole)CMBR anisotropies through the Sachs-Wolfe effect47.The observations imply48δρm1012GeV 7π|<0.51012GeV12.(13)Since−π<α1<+πis the a-priori range ofα1values and no particular value is preferred over any other,|α1π|=4·10−3,for example,f a may be as large as1016GeV.The presence of isocurvature perturbations constrains the smallα1scenario in two ways44. First,it makes it impossible to haveα1arbitrarily small sinceα1>δα1=H I1012GeV 5Second,one must require axion isocurvature perturbations to be consistent with CMBR bining Eqs.(11)and(12),and settingΩCDM=0.22,δρm4a f a24.(16) Let us keep in mind that the bounds(15)and(16)pertain only if T RH<T PQ.One may,for example,haveΩa=0.22,f a≃1012GeV,andΛI≃1016GeV,provided T RH>∼1012GeV,which is possible if reheating is sufficiently efficient.4Dark matter axion detectionAn electromagnetic cavity permeated by a strong static magneticfield can be used to detect galactic halo axions51.The relevant coupling is given in Eq.(4).Galactic halo axions have velocitiesβof order10−3and hence their energies E a=m a+1πgγm aMin(Q L,Q a)=0.510−26Watt V7Tesla 2C gγ1cm3m aB20V V d3xǫ| Eω|2(18) where B0( x)is the static magneticfield, Eω( x)e iωt is the oscillating electricfield andǫis the dielectric constant.Because the axion mass is only known in order of magnitude at best,the cavity must be tunable and a large range of frequencies must be explored seeking a signal.The cavity can be tuned by moving a dielectric rod or metal post inside it.For a cylindrical cavity and a homogeneous longitudinal magneticfield,C=0.69for the lowest TM mode.The form factors of the other modes are much smaller.The resonant frequency of the lowest TM mode of a cylindrical cavity is f=115MHz 1mKSVZ coupled axions at the nominal halo density of63ρa=7.510−25g/cm3over the mass range1.90<m a<3.35µeV60,61.The ADMX experiment is presently being upgraded toreplace the HEMT(high electron mobility transistors)receivers we have used so far with SQUID microwave amplifiers.HEMT receivers have noise temperature T n∼3K65whereas T n∼0.05K was achieved with SQUIDs66.In a second phase of the upgrade,the experiment willbe equipped with a dilution refrigerator to take full advantage of the lowered electronic noise temperature.When both phases of the upgrade are completed,the ADMX detector will have sufficient sensitivity to detect axions at even a fraction of the halo density.The ADMX experiment is equipped with a high resolution spectrometer which allows usto look for narrow peaks in the spectrum of microwave photons caused by discreteflows,or streams,of dark matter axions in our neighborhood.In many discussions of cold dark matter detection it is assumed that the distribution of CDM particles in galactic halos is isothermal. However,there are excellent reasons to believe that a large fraction of the local density of cold dark matter particles is in discreteflows with definite velocities64.Indeed,because CDM has very low primordial velocity dispersion and negligible interactions other than gravity,the CDM particles lie on a3-dim.hypersurface in6-dim.phase-space.This implies that the velocity spectrum of CDM particles at any physical location is discrete,i.e.,it is the sum of distinct flows each with its own density and velocity.We searched for the peaks in the spectrum of microwave photons from axion to photonconversion that such discreteflows would cause in the ADMX detector.We found none andplaced limits62on the density of any localflow of axions as a function of theflow velocitydispersion over the axion mass range1.98to2.17µeV.Our limit on the density of discreteflowsis approximately a factor three more severe than our limit on the total local axion dark matterdensity.AcknowledgmentsThis work is supported in part by the U.S.Department of Energy under grant DE-FG02-97ER41029.I am grateful to R.Woodard and A.Riotto for helpful discussions on topics related to axion isocurvature perturbations.References1.S.Adler,Phys.Rev.177(1969)2426;J.S.Bell and R.Jackiw,Nuov.Cim.60A(1969)47.2.S.Weinberg,Phys.Rev.D11(1975)3583.3.C.A.Baker et al.,hep-ex/0602020.4.M.Kobayashi and K.Maskawa,Progr.Theor.Phys.49(1973)652.5.R.D.Peccei and H.Quinn,Phys.Rev.Lett.38(1977)1440and Phys.Rev.D16(1977)1791.6.S.Weinberg,Phys.Rev.Lett.40(1978)223;F.Wilczek,Phys.Rev.Lett.40(1978)279.7.C.Vafa and E.Witten,Phys.Rev.Lett.53(1984)535.8.J.Ellis and M.K.Gaillard,Nucl.Phys.B150(1979)141.9.H.Georgi and L.Randall,Nucl.Phys.B276(1986)241.10.E.Witten,Phys.Lett.B149(1984)351.11.K.Choi and J.E.Kim,Phys.Lett.B154(1985)393;ibidem B165(1985)689.12.P.Svrcek and E.Witten,hep-th/0605206.13.J.Kim,Phys.Rev.Lett.43(1979)103;M.A.Shifman,A.I.Vainshtein and V.I.Zakharov,Nucl.Phys.B166(1980)493.14.M.Dine,W.Fischler and M.Srednicki,Phys.Lett.B104(1981)199;A.P.Zhitnitskii,Sov.J.Nucl.31(1980)260.15.J.E.Kim,Phys.Rev.D31(1985)1733.16.S.Weinberg in ref.6;W.A.Bardeen and S.-H.H.Tye,Phys.Lett.B74(1978)229;J.Ellis and M.K.Gaillard,Nucl.Phys.B150(1979)141;T.W.Donnelly et al.,Phys.Rev.D18(1978)1607;D.B.Kaplan,Nucl.Phys.B260(1985)215;M.Srednicki,Nucl.Phys.B260(1985)689;P.Sikivie,in’Cosmology and Particle Physics’,ed. E.Alvarez et al., World Scientific,1987,pp143-169.17.Axion reviews include:J.E.Kim,Phys.Rep.150(1987)1;H.-Y.Cheng,Phys.Rep.158(1988)1;R.D.Peccei,in’CP Violation’,ed.by C.Jarlskog,World Scientific Publ., 1989,pp503-551;M.S.Turner,Phys.Rep.197(1990)67;G.G.Raffelt,Phys.Rep.198 (1990)1.18.M.I.Vysotsky,Ya.B.Zeldovich,M.Yu.Khlopov and V.M.Chechetkin,JETP Lett.27(1978)502;D.Dicus,E.Kolb,V.Teplitz and R.Wagoner,Phys.Rev.D18(1978)1829and Phys.Rev.D22(1980)839.19.G.Raffelt and D.Dearborn,Phys.Rev.D36(1987)2211.20.D.Dearborn,D.Schramm and G.Steigman,Phys.Rev.Lett.56(1986)26.21.J.Ellis and K.Olive,Phys.Lett.B193(1987)525;G.Raffelt and D.Seckel,Phys.Rev.Lett.60(1988)1793;M.Turner,Phys.Rev.Lett.60(1988)1797;H.-T.Janka et al., Phys.Rev.Lett.76(1996)2621;W.Keil et al..Phys.Rev.D56(1997)2419.22.D.J.Gross,R.D.Pisarski and L.G.Yaffe,Rev.Mod.Phys.53(1981)43.23.L.Abbott and P.Sikivie,Phys.Lett.B120(1983)133;J.Preskill,M.Wise and F.Wilczek,Phys.Lett.B120(1983)127;M.Dine and W.Fischler,Phys.Lett.B120 (1983)137.24.R.Davis,Phys.Rev.D32(1985)3172and Phys.Lett.B18025.D.Harari and P.Sikivie,Phys.Lett.B195(1987)361.26.A.Vilenkin and A.E.Everett,Phys.Rev.Lett.48(1982)1867.27.P.Sikivie in Where are the elementary particles,Proc.of the14th Summer School onParticle Physics,Gif-sur-Yvette,1982,edited by P.Fayet et al.(Inst.Nat.Phys.Nucl.Phys.Particules,Paris,1983).28.P.Sikivie,Phys.Rev.Lett.48(1982)1156.29.S.Chang,C.Hagmann and P.Sikivie,Phys.Rev.D59(1999)023505.30.A.Vilenkin and T.Vachaspati,Phys.Rev.D35(1987)1138;R.L.Davis and E.P.S.Shellard,Nucl.Phys.B324(1989)167;A.Dabholkar and J.Quashnock,Nucl.Phys.B333(1990)815.31.C.Hagmann and P.Sikivie,Nucl.Phys.B363(1991)247.32.R.A.Battye and E.P.S.Shellard,Nucl.Phys.B423(1994)260,Phys.Rev.Lett.73(1994)2954and erratum-ibid.76(1996)2203.33.M.Yamaguchi,M.Kawasaki and J.Yokoyama,Phys.Rev.Lett.82(1999)4578.34.C.Hagmann,S.Chang and P.Sikivie,Phys.Rev.D63(2001)125018.35.D.Lyth,Phys.Lett.B275(1992)279;36.M.Nagaswa and M.Kawasaki,Phys.Rev.D50(1994)4821.37.A.Vilenkin,Phys.Rev.Lett.48(1982)59.38.P.J.Steinhardt and M.S.Turner,Phys.Lett.B129(1983)51;zarides,R.Schaefer,D.Seckel and Q.Shafi,Nucl.Phys.B346(1990)193.39.W.G.Unruh and R.M.Wald,Phys.Rev.D32(1985)831;M.S.Turner,Phys.Rev.D32(1985)843;T.DeGrand,T.W.Kephart and T.J.Weiler,Phys.Rev.D33(1986) 910;M.Hindmarsh,Phys.Rev.D45(1992)1130.40.D.B.Kaplan and K.M.Zurek,Phys.Rev.Lett.96(2006)041301.41.C.T.Hill and G.G.Ross,Nucl.Phys.B311(1988)253.42.C.J.Hogan and M.J.Rees,Phys.Lett.B205(1988)228;E.Kolb and achev,Phys.Rev.Lett.71(1993)3051,Phys.Rev.D49(1994)5040,and Ap.J.460(1996) L25.43.M.Axenides,R.H.Brandenberger and M.S.Turner,Phys.Lett.B126(1983)178;P.J.Steinhardt and M.S.Turner,Phys.Lett.B129(1983)51;A.D.Linde,Phys.Lett.B158 (1985)375;D.Seckel and M.Turner,Phys.Rev.D32(1985)3178;D.H.Lyth,Phys.Lett.B236(1990)408;44.M.Turner and F.Wiczek,Phys.Rev.Lett.66(1991)5.45.N.D.Birrell and P.C.W.Davies,Quantum Field Theory in Curved Space-Time,CambridgeUniversity Press,1982;A.Vilenkin and L.H.Ford,Phys.Rev.D25(1982)175;A.D.Linde,Phys.Lett.B108(1982)389;A.Starobinsky,Phys.Lett.B117(1982)175. 46.P.Crotty,J.Garcia-Bellido,J.Lesgourgues and A.Riazuelo,Phys.Rev.Lett.91(2003)17130147.R.Sachs and A.Wolfe,Ap.J.147(1967)73;P.J.E.Peebles,Ap.J.263(1982)L1;L.F.Abbott and M.B.Wise,Ap.J.282(1984)L47.48.S.Dodelson,Modern Cosmology,Academic Press,2003.49.S.-Y.Pi,Phys.Rev.Lett.52(1984)1725.50.D.N.Spergel et al.,astro-ph/0603449.51.P.Sikivie,Phys.Rev.Lett.51(1983)1415and Phys.Rev.D32(1985)2988.52.L.Krauss,J.Moody,F.Wilczek and D.Morris,Phys.Rev.Lett.55(1985)1797.53.C.Hagmann et al.,Rev.Sci.Inst.61(1990)1076.54.C.Hagmann,Ph.D.thesis,unpublished.55.S.DePanfilis et al.,Phys.Rev.Lett.59(1987)839and Phys.Rev.40(1989)3153.56.C.Hagmann et al.,Phys.Rev.D42(1990)1297.57.S.Matsuki and K.Yamamoto,Phys.Lett.B263(1991)523;S.Matsuki,I.Ogawa,K.Yamamoto,Phys.Lett.B336(1994)573;I.Ogawa,S.Matsuki and K.Yamamoto,Phys.Rev.D53(1996)R1740;A.Kitagawa,K.Yamamoto and S.Matsuki,hep-ph/9908445.58.C.Hagmann et al.,Phys.Rev.Lett.80(1998)2043.59.S.J.Asztalos et al.,Phys.Rev.D64(2001)092003.60.S.J.Asztalos et al.,Ap.J.Lett.571(2002)L27.61.S.J.Asztalos et al.,Phys.Rev.D69(2004)011101(R).62.L.Duffy et al.,Phys.Rev.Lett.95(2005)091304.63.E.I.Gates,G.Gyuk and M.S.Turner,Astroph.J.449(1995)123.64.J.R.Ipser and P.Sikivie,Phys.Lett.B291(1992)288;A.Natarajan and P.Sikivie.Phys.Rev.D72(2005)083513.65.R.F.Bradley,Nucl.Phys.B Proc.Suppl.72(1999)137.66.M.M¨u ck et al.,Appl.Phys.Lett.72(1998)2885;Nucl.Phys.B Proc.Suppl.72(1999)145.。

中药复方治疗急性胰腺炎作用机制的研究进展牛小龙1，2，姚广涛1，31 上海中医药大学研究生院，上海201203；2 上海中医健康服务协同创新中心；3 上海中医药大学创新中药研究院摘要：急性胰腺炎（AP）是临床常见的一种急腹症。

大多数AP患者为轻症，病程具有自限性，通常1～2周即可恢复。

但约20% AP患者会发展为重症急性胰腺炎（SAP），病死率为20%～40%。

西医治疗AP易引起继发性感染、腹膜炎、休克等并发症，整体治疗效果并不理想。

中医认为，AP起因于诸多病邪，包括热、湿、水、气、瘀等壅阻于胰、肝、胆、胃、脾、肠等脏腑，在治疗上应以“攻下通腑”“疏肝退热”“清热解毒”为突破点。

常用的中药复方包括大承气汤、大柴胡汤、大黄牡丹汤、柴芩承气汤、清胰汤等，其作用机制包括改善胃肠功能，修复肠黏膜屏障；抑制炎症反应，提高免疫功能；促进胰腺微循环；诱导胰腺腺泡细胞凋亡等。

这些中药复方以其多组分、多途径、多靶点相互作用，协同发挥治疗作用。

关键词：急性胰腺炎；中药复方；作用机制doi：10.3969/j.issn.1002-266X.2024.01.023中图分类号：R657.5+1 文献标志码：A 文章编号：1002-266X（2024）01-0093-05急性胰腺炎（AP）是临床常见的消化系统急症之一。

大多数AP患者为轻症，病程具有自限性，通常1～2周即可恢复。

但仍有约20% AP患者会发展为重症急性胰腺炎（SAP），病死率为20%～40%［1］。

西医治疗AP的主要方法包括立即禁食水、持续胃肠减压、静脉输液支持、抑制胃酸和胰液分泌等［2］。

但西医治疗易引起继发性感染、腹膜炎、休克等并发症，整体治疗效果并不理想。

中医药以其多组分、多途径、多靶点相互作用，协同发挥治疗作用，在治疗AP方面具有独特优势。

经典中药复方大承气汤、清胰汤能够减轻胰腺炎症，抑制病情加重［3］。

此外，大柴胡汤、大黄牡丹汤、柴芩承气汤等中药复方亦能通过改善胃肠功能、修复肠黏膜屏障、诱导细胞凋亡等基金项目：上海市科技计划项目资助（22S21901300）。