Monte Carlo simulation of the experiment MAMBO I and correction of neutron lifetime result
蒙特卡罗算法
xi 1 6xi mod11, ui+1 xi 1 /11 （a 6, m 11 ）
如果令 a 3, x0 1 ,得到序列: 1,3,9,5, 4,1,3,9........ 如果令 a 3, x0 2, 得到序列： 2, 6, 7,10,8, 2, 6.......
14
3. 指数分布 指数分布随机变量X的概率密度为
e x , x 0 f ( x) 0, x 0
指数分布常用来描述寿命问题.
15
二.随机数的生成
1.蒙特卡罗模拟的关键是生成优良的随机数。 2.在计算机实现中,我们是通过确定性的算法生成 随机数,所以这样生成的序列在本质上不是随机 的,只是很好的模仿了随机数的性质(如可以通过 统计检验)。我们通常称之为伪随机数(pseudorandom numbers)。 3.在模拟中，我们需要产生各种概率分布的随机数， 而大多数概率分布的随机数产生均基于均匀分布 U(0,1)的随机数。
2
基本思想很早以前就被人们所发现和 利用。17世纪，人们就知道用事件发 生的“频率”来决定事件的“概率”。 19世纪人们用投针试验的方法来决定π。 高速计算机的出现，使得用数学方法 在计算机上大量模拟这样的试验成为 可能。
3
从Buffon(蒲丰)投针问题谈起
Buffon 投针问题:平面上画很多平行线，间距为 a．向此平面投掷长为 l ( l < a) 的 针,此针与任一平行线相交的概率 p。
Monte Carlo Simulation Methods (蒙特卡罗模拟方法)
主要内容: 一. M-C方法概述. 二. 随机数的生成. 三. 模拟训练. 四. 实验题目.
成信院数学与信息科学系 李胜坤
MonteCarlo方法
2. 频率检验
检验每组观测频数 ni与理论频数mi ＝ N １/k之间相差的显著性
Monte Carlo 方法解决实际问题的过程中 ， 主要有以下几个内容 ① 建立简单而又便于实现的概率统计模型 ， 使所求的解正是该模型的某 一事件的概率或数学期望 ， 或该模型能够直接描述实际的随机过程。
② 根据概率统计模型的特点和计算的需求 ， 改进模型 ， 以便减小方差和 减低费用 ， 提高计算效率。
选择递推函数必须注意以下几点 ：
① 随机性好 ； ② 在计算机上容易实现 ； ③ 省时 ； ④ 伪随机数的周期长。
2 伪随机数的产生方法
最基本的伪随机数是均匀分布的伪随机数。 该方法是首先给一个 ２r位的数 ， 取其中间的 r位数码作为第一个伪随机数 ， 然后将这个数平方 ， 构成一个新的 ２ r位的数 ， 再取中间的 r位数作为第 二个伪随机数。 如此循环可得到一个伪随机数序列。 该方法的递推公式为
提高精度一位数 ， 抽样次数要增加１００ 倍 ； 减小随机变量的标准 差 ， 可以减小误差 。
Monte Carlo 方法具有以下四个重要特征 ：
① 由于 Monte Carlo 方法是通过大量简单的重复抽样来实现的 ， 因 此 ， 方法和程序的结构十分简单 。
② 收敛速度比较慢 ， 因此 ， 较适用于求解精度要求不高的问题。
具体试验步骤如下 ：
（1） 产生服从给定分布函数 g（x）的随机变量值 xi （2）检查 xi是否落入积分区域（a≤ x ≤ b） ， 如果满足条件 ， 则记录一次。
蒙特卡罗方法 (Monte Carlo simulation)教材
2019/4/15
Monte Carlo模拟
9
2.Monte Carlo方法简史
Stanislaw Ulam (1909-1984)
S. Ulam is credited as the inventor of Monte Carlo method in 1940s, which solves mathematical problems using statistical sampling.
1
Monte Carlo模拟
第一章 引言 （Introduction)
1. 2. 3. 4. Monte Monte Monte Monte Carlo方法 Carlo方法简史 Carlo模拟的应用 Carlo算法的主要组成部分
2019/4/15
Monte Carlo模拟
2
1.Monte Carlo方法
2019/4/15
Monte Carlo模拟
5
2.Monte Carlo方法简史 Buffon投针实验
1768年，法国数学家Comte de Buffon利用投针实验估 计的值
L
d
p
2019/4/15
d
2L
Monte Carlo模拟
6
2.Monte Carlo方法简史
Problem of Buffon’s needle: If a needle of length l is dropped at random on the middle of a horizontal surface ruled with parallel lines a distance d > l apart, what is the probability that the needle will cross one of the lines?
蒙特卡罗方法(MC)
蒙特卡罗方法（MC）蒙特卡罗（Monte Carlo）方法：蒙特卡罗（Monte Carlo）方法，又称随机抽样或统计试验方法，属于计算数学的一个分支，它是在本世纪四十年代中期为了适应当时原子能事业的发展而发展起来的。
传统的经验方法由于不能逼近真实的物理过程，很难得到满意的结果，而蒙特卡罗方法由于能够真实地模拟实际物理过程，故解决问题与实际非常符合，可以得到很圆满的结果。
这也是我们采用该方法的原因。
蒙特卡罗方法的基本原理及思想如下：当所要求解的问题是某种事件出现的概率，或者是某个随机变量的期望值时，它们可以通过某种“试验”的方法，得到这种事件出现的频率，或者这个随机变数的平均值，并用它们作为问题的解。
这就是蒙特卡罗方法的基本思想。
蒙特卡罗方法通过抓住事物运动的几何数量和几何特征，利用数学方法来加以模拟，即进行一种数字模拟实验。
它是以一个概率模型为基础，按照这个模型所描绘的过程，通过模拟实验的结果，作为问题的近似解。
可以把蒙特卡罗解题归结为三个主要步骤：构造或描述概率过程；实现从已知概率分布抽样；建立各种估计量。
蒙特卡罗解题三个主要步骤：构造或描述概率过程：对于本身就具有随机性质的问题，如粒子输运问题，主要是正确描述和模拟这个概率过程，对于本来不是随机性质的确定性问题，比如计算定积分，就必须事先构造一个人为的概率过程，它的某些参量正好是所要求问题的解。
即要将不具有随机性质的问题转化为随机性质的问题。
实现从已知概率分布抽样：构造了概率模型以后，由于各种概率模型都可以看作是由各种各样的概率分布构成的，因此产生已知概率分布的随机变量（或随机向量），就成为实现蒙特卡罗方法模拟实验的基本手段，这也是蒙特卡罗方法被称为随机抽样的原因。
最简单、最基本、最重要的一个概率分布是(0,1)上的均匀分布（或称矩形分布）。
随机数就是具有这种均匀分布的随机变量。
随机数序列就是具有这种分布的总体的一个简单子样，也就是一个具有这种分布的相互独立的随机变数序列。
蒙特卡罗方法 (Monte Carlo simulation)
[0, )
• The random vector is uniformly distributed on the region [0,d)×[0,). Accordingly, it has probability density function 1/d. • The probability that the needle will cross one of the lines is given by the integral
2014-8-27
Monte Carlo模拟
18
4.Monte Carlo算法的主要组成部分
Monte Carlo算法的主要组成部分 概率密度函数(pdf) 必须给出描述一个物理系统的一组概率密度函数;
随机数产生器 能够产生在区间[0,1]上均匀分布的随机数 抽样规则 如何从在区间[0,1]上均匀分布的随机数出发,随机抽 取服从给定的pdf的随机变量;
p
2014-8-27
0
l sin
0
1 d
2l dAd d
Monte Carlo模拟 8
2.Monte Carlo方法简史 Enrico Fermi
• 1930年,利用Monte Carlo方法研究中子的扩散 • 并设计了一个Monte Carlo机械装置，Fermiac,用于计算核 反应堆的临界状态
2014-8-27
Monte Carlo模拟
12
Monte Carlo模拟
第一章 引言 （Introduction)
1. 2. 3. 4. Monte Monte Monte Monte Carlo方法 Carlo方法简史 Carlo模拟的应用 Carlo算法的主要组成部分
蒙特卡洛模拟
蒙特卡洛模拟风险分析是我们制定的每个决策的一部分。
我们一直面对着不确定，不明确和变异。
甚至我们无法获得信息，我们不能准确的预测未来。
蒙特卡洛模拟( Monte Carlo simulation)让您看到了您决策的所有可能的输出，并评估风险，允许在不确定的情况下制定更好的决策。
什么是蒙特卡洛模拟( Monte Carlo simulation)蒙特卡洛模拟( Monte Carlo simulation)是一种计算机数学技术，允许人们在定量分析和决策制定过程中量化风险。
这项技术被专家们用于各种不同的领域，比如财经，项目管理，能源，生产，工程，研究和开发，保险，石油&天然气，物流和环境。
蒙特卡洛模拟( Monte Carlo simulation)提供给了决策制定者大范围的可能输出和任意行动选择将会发生的概率。
它显示了极端的可能性-最的输出，最保守的输出-以及对于中间路线决策的最可能的结果。
这项技术首先被从事原子弹工作的科学家使用；它被命名为蒙特卡洛，摩纳哥有名的娱乐旅游胜地。
它是在二战的时候被传入的，蒙特卡洛模拟( Monte Carlo simulation)现在已经被用于建模各种物理和概念系统。
蒙特卡洛模拟( Monte Carlo simulation)是如何工作的蒙特卡洛模拟( Monte Carlo simulation)通过构建可能结果的模型-通过替换任意存在固有不确定性的因子的一定范围的值(概率分布)-来执行风险分析。
它一次又一次的计算结果，每次使用一个从概率分布获得的不同随机数集。
根据不确定数和为他们制定的范围，蒙特卡洛模拟( Monte Carlo simulation)能够在它完成计算前调用成千上万次的重复计算。
蒙特卡洛模拟( Monte Carlo simulation)产生可能结果输出值的分布。
通过使用概率分布，变量能够拥有不同结果发生的不同概率。
概率分布是一种用来描述风险分析的变量中的不确定性的更加可行的方法。
Monte Carlo Simulations of the SU(2) Vacuum Structure
a r X i v:h e p -l a t /9312002v 1 1 D e c 19931Monte Carlo Simulations of the SU(2)Vacuum StructureA.R.Levia ∗aCenter for Theoretical Physics,Laboratory for Nuclear Science and Department of Physics,Massachusetts Institute of Technology,Cambridge,Massachusetts 02139U.S.A.and Department of Physics,Boston University,590Commonwealth Ave.,Boston,Massachusetts 02215U.S.A.Lattice Monte Carlo simulations are performed for the SU(2)Yang Mills gauge theory in the presence of an Abelian background with external sources to obtain information on the effective potential.The goal is to investigate the lowest Landau mode that,in the continuum one-loop effective potential,is the crucial mode for instability.It is shown that also in the lattice formulation this lowest Landau mode plays a very peculiar role,and it is important for the understanding of the vacuum properties.1.INTRODUCTIONTo understand the IR behavior of non-Abelian gauge theory a non-perturbative framework is necessary.Therefore,lattice formulation is par-ticularly useful.However,the comparison be-tween one-loop perturbative expansion and lat-tice regularization might give important informa-tion.From this comparison we can learn about the trustworthiness of perturbative expansion.Furthermore,the loop expansion can provide in-dicative information for the lattice quantities.Despite the importance of Yang-Mills theories,the complete solution of non-Abelian gauge the-ories has yet to be found.In order to gain a bet-ter understanding of these theories,a necessary first step is the study of their vacuum structures.Nevertheless,the vacuum structure of such theo-ries has yet to be understood even in the simplest case,i.e.SU(2)without matter.Moreover,the infrared properties of such theories must be stud-ied systematically if we want to have some clue on confinement.Many authors have studied one-loop effective potential for SU(2)and the possible consequences [1,2].In addition,several Monte Carlo simula-tions have already been done with a wild range of techniques in 3and 4Dimensions[2,3].2Hδb 3(xδµ2−yδµ1).(2)2With this choice,the well known Savvidy result for the one-loop effective potential is obtainedE(H)=H248π2 ln gH2−i6πν(4)whereνare the eigenvalues of the second deriva-tive of the action with respect to the gaugefields, and the summation is over all the eigenvalues. This expression is a natural consequence of the saddle point approximation:V(H)=V classic−¯hδηδη+O(¯h2)(5)whereΩis the volume factor.The eigenvalue of this problem can be found by realizing that there is the same symmetry as the one of the Lan-dau levels problem,therefore we have exactly the same class of solutions.That yield to the eigen-valuesν=k2+(2n+1+2S z)gH(6) where k=k20+k2z,n=0,1,2,...and S z=±1is the gluon polarization.For the lowest landau mode,i.e.n=0,S z=−1,we haveν<0that give the imaginary part of the effective potential.Note that with Abelian background it is possible to obtain all eigenvalues positive only by insertion of non-gauge invariant terms in the TTICE RESULTSUntil now we have argued that this lowest Lan-dau mode plays a very particular role on the one loop perturbative expansion.This motivated us to understand what happened to this mode on lattice where the task is non-perturbative.To extract information on this mode on lattice we noted that,due to thefiniteness of the lat-tice k z,is quantized as well,with k z=2mπ/L z, where m is an integer.Therefore,the lowest in-homogeneous z-mode,m=1,becomes stable for lattices the extent of which in the z-direction is smaller thanL criticz=2πgH.(7)This enables us to search for the critical sizeL criticz.The homogeneous mode,m=0,which is always unstable,is eliminated by imposing a delta condition in the path integral.We perform Monte Carlo simulations that gen-erate a backgroundfield H in the z direction in the presence of an Abelian source of strength j. We use a heat bath updating procedure with periodic boundary conditions and for the com-putational technique we address to reference[2]. To eliminate the homogeneous mode m=0we force the Polyakov line in the z direction to take a fixed value different from zero.Monte Carlo sim-ulations are made with a lattice volume L3∗L z, where L is the size of the x,y,t directions.The finite effect due to L is not so crucial.In fact,for the observable that we are interested,it is suffi-cient to use L=12−16.Simulations are made by varying L z from4to50.The expectation value of the plaquette in the 1-2plane P[F12]as a function ofβand j is mea-sured.We monitor the quantityX=P[F12(β,0)]−P[F12(β,j)]3interested,these other contributions are smoothly variable functions,and thus they will not affect the critical behavior.We evaluated X for different values of βand j performing,for large j 4500sweeps after dis-charging 500for thermalization.For smaller j we increase the number of sweeps until a significant amount of datawas collected.Our data [2]show that there is no sign of the unstable mode away from the critical βregion (β=2.1−2.5).The situation changes dramat-ically in the critical region where the instability appears as a decrease of the vacuum energy con-tribution to the plaquette in the z direction.This effect becomes more evident in the presence of strong sources.The fact that the instability dis-appears outside of the scaling window is a strong evidence that the instability is a distinguishing feature of the continuum rather than the strong coupling vacuum.In order to study the IR prop-erty of the continuum theory we should remain in the region where the instability is manifest.The disappearance of the instability as β→0may give a clue to understanding the difference between the physics of the continuum and the strong coupling lattice theory.We systematically analyzed the critical region of βfor several values of j .The values of L critic z(β,j )are obtained by interpolating the kink of X and taking the median value.L critic zwas found to be dependent on the source j and on βin this region according to the renormal-ization group dependence.It is clear from our data that L critic zbelongs to the confinement phase and that there is good agreement with the renor-malization group equation.Our data follow the same dependence of the deconfinament transition [6]as should be for a dimensional scale length.Hence the ratio between the L critic and the de-confinament scale parameter L dec is independent of β.The relevant physical quantities must be obtained in the limit of vanishing of the induced field Q (j )→0.From our data we obtain L critic zL critic z (m )=T decm .(10)The simulations show that these modes followwith good approximation eq.(10).This is strong evidence for the harmonicity of this modes.This is a quite surprising result because we are in a re-gion of strong non-perturbative effects.It is stim-ulating to think in terms of a parallel between the situation described above and the integer Quan-tum Hall Effect.This similarity is based on the dual picture (z ↔t ,E ↔H ,etc...),and on the plateau structures for each mode.AcknowledgementsIt is a pleasure to thank Janos Polonyi and Suzhou Huang for several useful discussions.REFERENCES1.G.K.Savvidy,Phys.Lett.71B (1977)133;S.G.Matinyan and G.K.Savvidy Nucl.Phys.B134(1978)539;N.K.Nielsen and P.Olesen,Nucl.Phys.B144(1978)376;H.B.Nielsen and M.Ninomiya Nucl.Phys.B156(1979)1;J.Ambjørn and P.Olesen,Nucl.Phys.B170(1980)60and 265;J.Ambjørn,B.Felsager and P.Olesen,Nucl.Phys.B175(1980)349;L.Maiani,G.Martinelli,G.C.Rossi and M.Testa,Nucl.Phys.B273(1986)275.2. A.R.Levi and J.Polonyi,preprint CTP 2161.3.J.Ambjørn et al,Phys.Lett.245B (1990)575;P.Cea and L.Cosmai,Phys.Lett.249B (1990)114,264B (1991)415;Phys.Rev.D43(1991)620;preprint BARI-TH 129/93;H.D.Trot-tier and R.M.Woloshyn,Phys.Rev.Lett.70(1993)2053.4.L.F.Abbott,Nucl.Phys.B185(1981)189,and references therein.5.S.Huang and A.R.Levi,preprint BU-HEP-93-22.6.J.Kuti,J.Polonyi and K.Szlachanyi,Phys.Lett.98B (1981)199;E.Kovacs,Phys.Lett.118B (1982)125.。
Monte-Carlo模拟
曼哈顿计划 Buffon投针实验 大数定律
基本思想：当所求问题是某种随机事件出现的 概率，或者是某个随机变量的期望值时，通过 某种“试验”的方法，以这种事件出现的频率 估计该随机事件的概率，或者得到这个随机变 量的某些数字特征，并将其作为问题的解。
3) 建立各种估计量：一般说来，构造了 概率模型并能从中抽样（即实现模拟实 验）后，我们就要确定一个随机变量， 作为所求的问题的解。即针对模拟实验 的结果考察其统计特性（样本均值、方 差、置信区间等），建立各种估计量， 从中得到问题的解。
明确问题，建立模型收集和整理数据资料 编制程序，模拟运行分析模拟输出结果
逆变换法的具体步骤：
•确定随机变量X的概率分布函数F(x)；
例：指数分布的Biblioteka 布函数为：1 e x , x 0, F ( x) x 0. 0,
F ( x) 1 e , x 0
解得
可取
1 x ln(1 ), 1 x ln .
模拟的优点：简单、快速、适应性强
能相对容易地近似很复杂的随机系统，问题 的几何形状的复杂性对其影响不大；
可以在广泛的条件下估计候选方案的性能； 建模者可以在不同层次的水平上进行控制；
模拟的缺点： 建立和运行模拟模型可能相当昂贵； 模拟模型的随机性使得结论受到限制
2、随机数和随机变量的生成 2.1 均匀分布随机数的生成
n=input('输入模拟次数：'); count=0; for i=1:n, rt1=rand; %模拟随机变量t1（火车从A站出发的时刻） if rt1<0.7 T1=0; elseif rt1>=0.7 & rt1<0.9 T1=5; else T1=10; end T2=30+randn*2; %模拟随机变量t2（火车的运行时间） %模拟随机变量t3（他到达B站的时刻） rt3=rand; if rt3<0.3 T3=28; elseif rt3>=0.3 & rt3<0.7 T3=30; elseif rt3>=0.7 & rt3<0.9 T3=32; else T3=34; end if T3 < T1 + T2, count=count+1; end end%for prob=count/n
蒙特卡罗(Monte Carlo)方法简介
蒙特卡罗(Monte Carlo)方法简介蒙特卡罗(Monte Carlo)方法简介蒙特卡罗(Monte Carlo)方法，也称为计算机随机模拟方法，是一种基于"随机数"的计算方法。
一起源这一方法源于美国在第二次世界大战进研制原子弹的"曼哈顿计划"。
Monte Carlo方法创始人主要是这四位：Stanislaw Marcin Ulam, Enrico Fermi, John von Neumann（学计算机的肯定都认识这个牛人吧）和Nicholas Metropolis。
Stanislaw Marcin Ulam是波兰裔美籍数学家，早年是研究拓扑的，后因参与曼哈顿工程，兴趣遂转向应用数学，他首先提出用Monte Carlo方法解决计算数学中的一些问题，然后又将其应用到解决链式反应的理论中去，可以说是MC方法的奠基人；Enrico Fermi是个物理大牛，理论和实验同时都是大牛，这在物理界很少见，在“物理大牛的八卦”那篇文章里提到这个人很多次，对于这么牛的人只能是英年早逝了（别说我嘴损啊，上帝都嫉妒！）；John von Neumann可以说是计算机界的牛顿吧，太牛了，结果和Fermi一样，被上帝嫉妒了；Nicholas Metropolis，希腊裔美籍数学家，物理学家，计算机科学家，这个人对Monte Carlo方法做的贡献相当大，正式由于他提出的一种什么算法（名字忘了），才使得Monte Carlo方法能够得到如此广泛的应用，这人现在还活着，与前几位牛人不同，Metropolis很专一，他一生主要的贡献就是Monte Carlo方法。
蒙特卡罗方法的名字来源于摩纳哥的一个城市蒙地卡罗，该城市以赌博业闻名，而蒙特•罗方法正是以概率为基础的方法。
与它对应的是确定性算法。
二解决问题的基本思路Monte Carlo方法的基本思想很早以前就被人们所发现和利用。
早在17世纪，人们就知道用事件发生的"频率"来决定事件的"概率"。
蒙特卡洛方法2
由大数定律可知，当 n ，样本的均值趋 n 逐渐增大 向与理论分布的期望，因此利用样本容量 n 这一趋势来模拟 这一趋势，在这种趋势下， 应该呈现出越来 样本的均值与理论分布期望的误差 越小的趋势
一些人进行了实验，其结果列于下表 ：
实验者 沃尔弗(Wolf) 年份 1850 投计次数 5000 π的实验值 3.1596
斯密思(Smith)
福克斯(Fox) 拉查里尼 (Lazzarini)
1855
1894 1901
3204
1120 3408
3.1553
3.1419 3.1415929
例2. 射击问题（打靶游戏）
MANIAC the Computer and the Man
Seated is Nick Metropolis, the background is the MANIAC vacuum tube computer
Monte Carlo
Direct Monte Carlo Important sampling Monte Carlo Metropolis Monte Carlo Dynamic Monte Carlo
因此，可以通俗地说，蒙特卡罗方法是用随机试 验的方法计算积分，即将所要计算的积分看作服从某 种分布密度函数f(r)的随机变量ｇ(r)的数学期望
g g (r ) f (r )dr
0
通过某种试验，得到Ｎ个观察值r1，r2，…，rN（用概 率语言来说，从分布密度函数 f(r) 中抽取 Ｎ 个子样 r1 ， r2 ， … ， rN ，），将相应的 Ｎ 个随机变量的值 g(r1) ， g(r2)，…，g(rN)的算术平均值
蒙特卡罗(Monte Carlo)方法
/share/detail/5568877蒙特卡罗(Monte Carlo)方法，或称计算机随机模拟方法，是一种基于“随机数”的计算方法。
这一方法源于美国在第一次世界大战进研制原子弹的“曼哈顿计划”。
该计划的主持人之一、数学家冯·诺伊曼用驰名世界的赌城—摩纳哥的Monte Carlo—来命名这种方法，为它蒙上了一层神秘色彩。
Monte Carlo方法的基本思想很早以前就被人们所发现和利用。
早在17世纪，人们就知道用事件发生的“频率”来决定事件的“概率”。
19世纪人们用投针试验的方法来决定圆周率π。
本世纪40年代电子计算机的出现，特别是近年来高速电子计算机的出现，使得用数学方法在计算机上大量、快速地模拟这样的试验成为可能。
考虑平面上的一个边长为1的正方形及其内部的一个形状不规则的“图形”，如何求出这个“图形”的面积呢？Monte Carlo方法是这样一种“随机化”的方法：向该正方形“随机地”投掷N个点落于“图形”内，则该“图形”的面积近似为M/N。
可用民意测验来作一个不严格的比喻。
民意测验的人不是征询每一个登记选民的意见，而是通过对选民进行小规模的抽样调查来确定可能的优胜者。
其基本思想是一样的。
科技计算中的问题比这要复杂得多。
比如金融衍生产品（期权、期货、掉期等）的定价及交易风险估算，问题的维数（即变量的个数）可能高达数百甚至数千。
对这类问题，难度随维数的增加呈指数增长，这就是所谓的“维数的灾难”(Course Dimensionality)，传统的数值方法难以对付（即使使用速度最快的计算机）。
Monte Carlo方法能很好地用来对付维数的灾难，因为该方法的计算复杂性不再依赖于维数。
以前那些本来是无法计算的问题现在也能够计算量。
为提高方法的效率，科学家们提出了许多所谓的“方差缩减”技巧。
另一类形式与Monte Carlo方法相似，但理论基础不同的方法—“拟蒙特卡罗方法”(Quasi-Monte Carlo方法)—近年来也获得迅速发展。
数学建模专题三 Monte Carlo模拟
2019/5/18
Lxy, China Jiliang Universty
16
结果比较
理论计算和模拟结果的比较
分类 项目
无效射击
有效射击
数学建模专题三 -Monte Carlo模拟
平均值
模拟 理论
０．６５ ０．７５
０．３５ ０．２５
０．５ ０．３３
虽然模拟结果与理论计算不完全一致，但它却能更加真实地表 达实际战斗动态过程．
x=randperm(6); y=x(1); switch y
Lxy, China Jiliang Universty
10
问题分析
数学建模专题三 -Monte Carlo模拟
需要模拟出以下两件事：
[1] 观察所对目标的指示正确与否 模拟试验有两种结果，每一种结果出现的概率都是1/2． 因此，可用投掷一枚硬币的方式予以确定，当硬币出现正面时为
指示正确，反之为不正确．
4
举例
数学建模专题三 -Monte Carlo模拟
Buffon投针实验
1768年，法国数学家Comte de Buffon利用投针实验估计的值
2019/5/18
L
d
p
2L d
Lxy, China Jiliang Universty
5
Solution
数学建模专题三 -Monte Carlo模拟
... ...
2019/5/18
Lxy, China Jiliang Universty
19
Matlab中的取整函数
数学建模专题三 -Monte Carlo模拟
fix(x) : 截尾取整，直接将小数部分舍去 floor(x) : 不超过 x 的最大整数 ceil(x) : 不小于 x 的最小整数 round(x) : 四舍五入取整
高分子物理的Monte Carlo模拟
Distribution(%)
0.2 BLOCK A
0.4 BLOCK A 0.6 BLOCK A 0.8 BLOCK A
Size of Cores
20% 二嵌段共聚物的体系核尺 寸的分布
3.3、三嵌段共聚物的溶液 ABA 和BAB共聚物 其中 AB =1, AS =1,
BS=-1
10%
2. 统计物理学中的Monte Carlo 模拟方法，
K. Binder, D.W. Heermann著, 秦克诚译， 北京大学出版社，1994年.
0.2(6个链节)硬段
0.3(9个链节)硬段
0.4(12个链节)硬段
0.2硬段
0.3硬段
0.4硬段
0.2(6个链节)硬段
0.3(9个链节)硬段
0.4(12个链节)硬段
0.2硬段
0.3硬段
0.4硬段
参考文献
1. 杨玉良，张红东，高分子科学中的Monte Carlo 方法，复旦大学出版社，1993年.
第四部分：高分子物理的Monte Carlo 模拟
及键涨落模型
一 统计物理学中的Monte Carlo 模拟 1. 平衡态
任一物理量A的统计平均 值<A>
1 A Ae Z
H kT
dpdv
Z e
H kT
dpdv
Z为配分函数
哈密顿量 H
2 P H rU m2
U ( r ) ij [( ) ( ) ] rij rij
2 0 0 10 20 30 40
End-To-End Distance
22 20 18 16
End-To-End Distance
Monte Carlo Simulations of the Random-Field Ising Model
a r X i v :c o n d -m a t /0009228v 1 [c o n d -m a t .d i s -n n ] 14 S e p 2000Monte Carlo Simulations of the Random-Field Ising ModelW.C.Barber and D.P.BelangerDepartment of Physics,University of California,Santa Cruz,CA 95064,USAIsing Monte Carlo simulations of the random-field Ising system Fe 0.80Zn 0.20F 2are presented for H =10T.The specific heat critical behavior is consistent with α≈0and the staggered magnetization with β≈0.25±0.03.It has recently been shown experimentally [1]and through Monte Carlo (MC)simulations [2]that the three dimensional (d =3)dilute,anisotropic antiferromagnet Fe x Zn 1−x F 2in an applied uniform field exhibits the equilibrium critical behavior of the random-field Ising model (RFIM)[3]only for magnetic concentrations x >0.75.Although there is some agreement between previous MC and equilibrium experimental data,particularly the staggered susceptibility critical exponent,where γ=1.7±0.2for MC [4]and γ=1.58±0.08for experiment [1],and the correlation length critical exponent,where ν=1.1±0.2for MC [4]and ν=0.88±0.05for experiment [1],the experimental value [5]of the specific heat exponent,α=0.00±0.02,disagrees substantially with the simulation result [4]α=−0.5±0.2.The order parameter exponent,β=0.00±0.05from simulations [4],has yet to be measured experimentally.The MC exponents violate the Rushbrooke scaling inequality 2β+γ+α≥2,indicating possible inconsistencies.The previous MC studies employed finite scaling on ferromagnetic lattices of sizes L ×L ×L ,with L ≤16.Simulations were made over many random-field configurations and equilibrium was ensured by insisting that different starting configurations yielded consistent results at each temperature.In the present MC study,we have explored a much larger antiferromagnetic lattice and modeled our simulations closely after Fe 1−x Zn x F 2and the thermal cycling procedures used in experiments.Our large lattice size prevented the averaging over many random-field configurations,though the results obtained for a few different configurations yielded consistent results.We determined the staggered magnetization,M s ,and C m versus T .The resulting C m more closely mimics the symmetric,nearly logarithmic C m peak of the experiments [5]than the nondivergent cusp found in the earlier MC studies.In addition,we find β=0.25±0.03,in disagreement with the previous MC studies.The body-centered-tetragonal magnetic lattice of Fe 0.80Zn 0.20F 2is modeled as two cubic sub-lattices delineated as one-dimensional arrays bit coded to allow for a very large lattice size,2L ×L ×L with L =128,corresponding to 3.4×106spins.Periodic boundary conditions are imposed.The dilute MC Ising Hamiltonian is H =J 2 <ij>ǫi ǫj S i S j −h i ǫi S i (1)where S i =±2and ǫi =1if site i is occupied and zero if not,J 2=3.17K and h =8.73K.This approximates the experimental system which has one dominant Heisenberg exchange interaction,J 2,and a large single-ion anisotropy.The very small,frustrating interactions,J 1and J 3,are not included.The MC value for J 2is 0.60times that of FeF 2so that T c (0)corresponds roughly tothat of the experimental system.The value of h is chosen to have roughly the same effect as a field H =10T in the real system,as determined by the fields at which low temperature spin flips occur.This field also yields a shift T c (0)−T c (H )consistent with H =10T in the real system [6].Henceforth we refer to the applied uniform field as H =10T.For H =0,the sample is cooled to low temperatures in steps of 0.01K while magnetic sites are randomly visited an average of N times,where 2000<N <5000,and flipped with a probability given by the metropolis algorithm.Similarly,the sample is heated from an ordered state to above the transition.No hysteresis is observed at H =0using these two procedures.For H >0,two thermal cycling procedures are used to mimic those of the experiments.Each site is visited 5000<N <10000times per 0.01K temperature step.For zero-field cooling (ZFC),the sample is1started at H=0in a well ordered state(in experiments this state is prepared by cooling in H=0), thefield is subsequently applied,and the sample is heated through T c(H).Infield-cooling(FC),the sample is cooled with H applied.We observe no significant hysteresis with these two procedures, consistent with the experiments.M s and C m are calculated for each temperature.The simulations are performed on Linux/Gnu computers with CPU speeds between350to450MHz.Individual runs with L=128at5000MC steps per spin and covering10K in T typically take one week of CPU time.We have not seen any significant dependence on N in the range studied.Figure1shows C m versus T for L=128and N=5000for H=0and for ZFC and FC at H=10T.Each point represents a point-by-point derivative of the energy which is itself averaged over0.1K intervals.The behavior is very similar to that observed in birefringence and pulsed heat experiments[5].Namely,the H=0behavior exhibits the asymmetric cusp(α<0)characteristic of the random-exchange Ising model and,for H>0,a much more symmetric peak near T c(H) consistent with the symmetric,nearly logarithmic divergence.Previous MC studies employingfinite size scaling did not yield this behavior.Figure2shows the ZFC M s versus T−T c for L=128and N=5000for H=0and H=10T. Fits of the data for various ranges of|t|within10−3<|t|<10−1,yield the critical order-parameter exponentβ=0.35±0.04for H=0(REIM)andβ=0.25±0.03at H=10T(RFIM).The T c(H) values used in thefits were constrained by thefits to C m(H)To demonstrate that the results of the simulations are close to equilibrium,we did a simulation just below the transition for L=128,H=10T,and T=60.5.One lattice was started fully ordered,another fully ordered but with the sublattices reversed,and one randomly ordered.All three lattices converged to the same energy and M s within expectedfluctuations after5000steps per spin.Usingβ=0.25±0.03from MC andα=0.02±0.02andγ=1.58±0.08from experiments, we obtain2β+γ+α=2.08±0.14,which is consistent with the Rushbrooke scaling relation as an equality.Further studies will be needed to understand why the MC simulations with large an-tiferromagnetic lattices and thermal cycling yield different results than the earlier ferromagnetic MC simulations withfinite scaling analyses[4].For simulations of RFIM critical behavior employ-ing dilute antiferromagnets,the magnetic concentration should be kept well above the percolation threshold concentration for vacancies.This work has been supported by the Department of Energy Grant No.DE-FG03-87ER45324.FIG.1.C m vs.T for H=0and for ZFC(triangles)and FC(inverted triangles)at10T with L=128 and x=0.80.Note the symmetric peak for H=10T.The curves arefits to the data withα=−0.09and 0.00for the H=0and H=10T cases,respectively.FIG.2.M s vs.T−T c for H=0and ZFC at H=10T with L=128and x=0.80.The curves arefits to the data withβ=0.35and0.25for the H=0and H=10T cases,respectively.3。
Molecular Simulation5Monte Carlo Simulation
• Histogram the occupancy number for each state
– n1 = 3 – n2 = 5 – n3 = 4
p1 = 0.25 p2 = 0.42 p3 = 0.33
123
• After large enough steps, a limiting
p1(2n)
p
(n) 22
p
(n) 32
p1(3n)
p p
(n) 23
(n) 33
– probabilities of going from state i to j in exactly n step
defines
p
(n) ij
9
• Define pi(0) as a unit state vector
7
Distribution of State Occupancies
• Consider process of repeatedly moving from one state to the next based on P, choosing each subsequent state according to P
For NVT ensemble, the limiting distribution is already known: (e-βE1/Q, e –βE2/Q, e –βE3/Q, …………)
How to design the transition matrix Π?
p pP a p11a p 21b p31c b p 21a p 22b p32c c p31a p32b p33c
计算机模拟-17
the previous state 9. determine the value of the desired physical
quantities, e.g. spin average etc
1. initial state: (a) random distribution of spinup or spin-down; (b) alignment up or down
2. chose a site: (a) randomly or (b) sequentially 3. flip the spin at the site, (single spin flip) 4. calculate the energy Etrial, and compute E=
write random number generation routine takes times, and normally does not work properly
copy routine from publications; download it from computer web-sites
8
Application in Statistical Physics
An observable A(x)
Adx(xA )expH(N(x)/kBT) dexxpH(N(x)/kBT)
macroscopic system N 1022 atoms; Monte Carlo accessible N <103 - 106 atoms, sampling
计算机模拟在铁电 物理研究中的应用
Monte_Carlo算法模拟
一些人进行了实验，其结果列于下表 ：
实验者 沃尔弗(Wolf) 年份 1850 投计次数 5000 π的实验值 3.1596
斯密思(Smith)
福克斯(Fox) 拉查里尼 (Lazzarini)
1855
1894 1901
3204
1120 3408
3.1553
3.1419 3.1415929
蒙特卡罗投点法是蒲丰投针实验的推广： 在一个边长为a的正方形内随机投点， 该点落在此正方形的内切圆中的概率 y (a/2,a/2) 应为该内切圆与正方形的面积比值， 即 π a/2 2 : a2 π/4 n=10000; a=2; m=0; o x for i=1:n x=rand(1)*a; y=rand(1)*a; if ( (x-a/2)^2+(y-a/2)^2 <= (a/2)^2 ) m=m+1; end end disp(['投点法近似计算的π为: ',num2str(4*m/n)]);
rand('seed',0.1); rand(1) ％每次运行程序产生的值是相同的 rand('state',sum(100*clock)*rand);
rand(1) %每次重新启动matlab时，输出的随机数不一样
注意: 产生一个参数为λ的指数分布的随机数应输入 exprnd(1/λ)
产生m×n阶参数为A1,A2,A3的指定分布'name'的随机数矩阵 random('name',A1,A2,A3,m,n) 举例: 产生2×4阶的均值为0方差为1的正态分布的随机数矩阵 random('Normal',0,1,2,4) 'name'的取值可以是(详情参见help random)： 'norm' or 'Normal' / 'unif' or 'Uniform' 'poiss' or 'Poisson' / 'beta' or 'Beta' 'exp' or 'Exponential' / 'gam' or 'Gamma' 'geo' or 'Geometric' / 'unid' or 'Discrete Uniform' ……
薄膜生长的理论模型与MonteCarlo模拟
薄膜生长的理论模型与Monte C arlo 模拟Ξ杨 宁1) 陈光华1) 张 阳2) 公维宾1) 朱鹤孙2)1)(北京工业大学材料科学与工程学院,北京 100022)2)(北京理工大学材料科学研究中心,北京 100081)(2000年4月21日收到;2000年6月11日收到修改稿) 用Monte Carlo 方法模拟了薄膜的二维生长.引入Morse 作用势描述粒子间的相互作用情况,研究了相互作用范围α对薄膜生长初期形貌的影响.薄膜的生长经历了由临界核的形成、团的长大、形成迷津结构到连续成膜4个阶段.模拟结果表明,随着沉积粒子数的增加,成团粒子所占百分比下降.这与实际情况相符.Ξ国家自然科学基金(批准号:69876003和19874007)及北京市自然科学基金(批准号:2982013)资助的课题.关键词:薄膜生长,Monte Carlo 方法,Morse 势PACC :6855,31201 引言随着科学技术的发展,各行业对新材料的需求日益迫切.薄膜材料,特别是电子材料在现代材料中占有日益重要的地位.为了提高薄膜质量,人们总是通过调整制膜工艺条件来控制成膜的微观过程.因此要想尽快地探索出成膜的最佳工艺条件,就需要对粒子在衬底上的成膜过程和各种参数对成膜的影响有较深入的了解,而计算机模拟理论恰好提供了这样一个有效的手段.通过对不同模型的模拟及不同参数的设置,可以了解不同工艺条件下的成膜过程及各参数变化对成膜的影响.本文用Monte Carlo 方法模拟了二维平面上粒子的成膜情况.所考虑的成膜过程是大量微观粒子在给定宏观约束条件下的集体行为,但就每一个粒子的行为而言,则是随机的.结合粒子的随机运动,根据过程的物理特性设定一些概率规则,理论模拟结果与实际情况符合较好.2 薄膜生长的理论模型利用周期性边界条件建立一个100×100的二维方格点阵.考虑到衬底点阵原子势能的极小值应当位于各阵点的正上方,所以落下的粒子只能位于某个阵点所在的坐标处,而不会落在阵点之间的位置.这相当于按衬底的晶格结构外延生长,仍生成二维结构的膜.这样,以点阵常数为单位,粒子坐标只取整数值.同理,粒子扩散的步长也为点阵常数.粒子之间的作用势采用讨论材料物理性质时常用的Morse 势,即<(r ij )=D{exp [-2α(r ij /r 0-1)]-2exp [-α(r ij /r 0-1)]},(1)r ij 为编号为i ,j 两粒子之间的距离;D 和α分别为相互作用的强度和范围;r 0为势能最低点的位置,即平衡位置;D 可用下列参数表示:E ff ,E fs ,E ss ;E ff 为成膜粒子间的结合能;E fs 为成膜粒子与衬底粒子间的结合能;E ss 为衬底粒子间的结合能.本文只考虑在二维情况下的单层膜的生长.粒子在衬底上的运动可分为三种情况:(1)吸附:粒子被吸附在衬底表面.在本模型中假设为物理吸附;(2)扩散:粒子在衬底上迁徙,粒子的扩散方向应使系统的整体能量降低;(3)蒸发:粒子获得一足够大的能量使其脱离衬底.本文只涉及前两种情况,粒子的再蒸发将另文讨论.在本模型中,每次有一个粒子随机落到衬底上,若粒子落到第二层上则根据下落的位置,按照一定概率扩散到衬底上.用(1)式计算该点处粒子的势能,在该粒子的4个最近邻点(图1中的1,2,3,4格点)中,选择可以使系统能量处于更低的点扩散,若无这样的点则根据Boltzmann 因子ρB 判断是否扩散:第49卷第11期2000年11月100023290/2000/49(11)/2225205物 理 学 报ACTA PHYSICA SIN ICAVol.49,No.11,November ,2000ν2000Chin.Phys.S oc.ρB =e-ΔH/k T,(2)ΔH 为系统能量的改变量;k 为Boltzmann 常数;T为系统温度.因此当新位置的能量高于原来位置的能量时,新位置就有一定的概率被接受.这样处理可以保证当系统处于一个局部范围内的能量最低点时,系统有可能离开该极小值点,到达整个系统能量的最小值处,否则系统将陷于该极小值点处,无法使系统的自由能达到最小.图1 沉积粒子及其4个最近邻位置本文采用了周期性边界条件,这样当粒子从衬底左边扩散出去后,相应地,从右边相对位置会进来一个粒子;当粒子从上边出去后,会从下边相应位置进来一个粒子.同理,当左右、上下相对位置分别有粒子存在时,认为这些粒子相应的近邻处有粒子存在,统计成团个数时将其计算在内. 计算某格点上粒子的势能时,并不只是考虑4个最近邻处粒子对其的影响,而是根据Morse 势中作用范围α的取值来决定.如图2所示,根据不同的α值,可取得不同的作用范围R ,在R 范围内的粒子都被予以考虑.图2 Morse 势势能曲线根据上述模型编写了程序,实时模拟了薄膜的二维生长.3 结果与讨论计算了在方形格子表面上薄膜的沉积过程,表面格点数为L 2=100×100.图3(a )—(f )为参数为α=2,T =300K ,E ff =012eV ,E fs =014eV 时薄膜的形貌图.可以看出在此条件下,粒子趋向于成团生长,生长过程分别经历了临界核的形成、原子团的长大、迷津结构的形成和连续成膜[11]等几个阶段.图4(a )—(f )为其他条件相同,而α=6时薄膜的形貌图.可以看出与α=2时有很大区别.此时沉积的粒子分布十分分散,不易成团生长.(a )N =100(b )N =400(c )N =700(d )N =1000(e )N =4000图3 α=2时薄膜形貌图(f )N =70006222物 理 学 报49卷(a)N=100(b)N=400(c)N=700 (d)N=1000(e)N=4000图4 α=6时薄膜形貌图(f)N=7000 图5为成团个数C随沉积粒子数N变化的规律.可以看出,团的个数随沉积粒子数的增加而迅速增加,到达极大值后又逐步下降.对应于α=2时,成团数目增加迅速,极值点出现在沉积粒子数为1000—2000之间,且在下降阶段的末期趋近于1.对于α=6时,团的数目开始时很少,经过阈值后,团的数目增加较快,但与α=2时相比,仍较为缓慢.α=6时曲线的极大值点出现在3000—4000之间.这参照图3、图4可以很容易地得到解释.图5中的极大值表示稳定的原子团密度达到了饱和,而过了极大值,由于原子团之间的合并占主要地位,使原子团的数目减小.趋近于1的转折点,对应于原子团合并成一个大团,形成迷津结构,这时沉积到衬底上的粒子主要是填补迷津结构的空隙,这标志着连续成膜的开始. 图6为成团粒子数占沉积粒子数的百分比ρc随沉积粒子数N变化的情况.可以看到α=2时,ρc随沉积粒子数的增加而急剧上升,随后又缓慢下降,极值点约在1000处.α=6时,上升趋势则相对缓慢得多,极值点出现在5000左右.值得注意的是,两条曲线在到达极值点后都缓慢下降,而并不是趋近于饱和[2].这表明在此阶段沉积的粒子中用于在衬底上成团的粒子数所占的比例减小了.这是由于对模型的不同处理方法造成的.在本模型中,对于落到第二层上的粒子,若其落到团的边缘,则可扩散到下一层.若该粒子的4个可扩散方向上都有粒子存在时,则认为它所落的位置靠近团的中央,此时扩散到下一层的概率极小,因此,认为它不再扩散而是落在第二层上.在总的沉积粒子数中记入该点,但在计算第一层的成核情况时,并不记入该点.由此可知曲线的下降趋势说明第二层以上粒子的沉积所占的比重逐渐增加,这在三维生长模型中得到了证实.在其他的二维模型[2—4]中,若有落到第二层上的粒子,则取消重落或让其在第二层上扩散直至跃迁到第一层,在这种情况下,所有沉积的粒子均用于第一层上的成膜,因此最后会趋近饱和.722211期杨 宁等:薄膜生长的理论模型与Monte Carlo模拟图5 成团个数随沉积粒子数变化规律图6 成团粒子数占沉积粒子数的百分比ρc 随沉积粒子数N 变化规律4 结论11用计算机模拟研究了薄膜生长机理,给出了微观生长过程的详细图像.21Morse 势中相互作用范围α影响薄膜的形貌,α=2时,薄膜趋向于成团生长,α=6时,薄膜趋向于分散生长.31得到的ρc 2N 曲线有一下降趋势,此趋势说明了随着沉积粒子的增加多层生长所占比重逐渐增加,这更符合实际情况.[1]Z.Q.Xue ,Q.D.Wu ,H.Li ,Thin Film Physics ,1st ed.(Pub 2lishing House of Electronics Industry ,Beijing ,1991),p.20(in Chinese )[薛增泉、吴全德、李 浩,薄膜物理,第1版(电子工业出版社,北京,1991),第20页].[2]G.Yu ,G.Z.Yue ,J.Qi ngdao Instit ute of Chemical Technolo 2gy ,17(1996),279(in Chinese )[于 工、岳国珍,青岛化工学院学报,17(1996),279].[3]R.Smith ,A.Richter ,Thi n Soli d Fil ms ,3432344(1999),1.[4]Z.L.Liu ,H.L.Wei ,H.W.Wang ,J.Z.Wang ,Acta PhysicaSi nica ,48(1999),1302(in Chinese )[刘祖黎、魏合林、王汉文、王均震,物理学报,48(1999),1302].8222物 理 学 报49卷MONTE CAR LO SIMU LATION OF THIN FILMS GROWTH ΞY AN G N IN G 1) C HEN G UAN G 2HUA 1) Z HAN G Y AN G 2) G ON G W EI 2BIN 1) Z HU H E 2SUN 2)1)(School of M aterials Science and Engineering ,Beijing Polytechnic U niversity ,Beijing 100022,China )2)(Research Center of M aterials Science ,Beijing Instit ute of Technology ,Beijing 100081,China )(Received 21April 2000;revised manuscript received 11J une 2000)A BSTRACTTwo 2dimensional growth of thin films is simulated using the Monte Carlo method.The effect of the interaction be 2tween particles on the morphology of the films at the preliminary stage has been investigated using Morse potential.The growth of films consists of four stages ,including formation of critical nucleus ,growth of clusters ,formation of maze struc 2ture and appearance of continuous film.The simulated results indicate that the percentage of particles ,with monomers be 2ing not included ,decreases with the increase of the number of adatoms ,which is in agreement with the actual conditions.K eyw ords :thin film growth ,Monte Carlo method ,Morse potential PACC :6855,3120ΞProject supported by the National Natural Science Foundation of China (Grant Nos.69876003and 19874007)and the Natural Science Founda 2tion of Beijing ,China (Grant No.2982013).922211期杨 宁等:薄膜生长的理论模型与Monte Carlo 模拟。
- 1、下载文档前请自行甄别文档内容的完整性，平台不提供额外的编辑、内容补充、找答案等附加服务。
- 2、"仅部分预览"的文档,不可在线预览部分如存在完整性等问题,可反馈申请退款(可完整预览的文档不适用该条件!)。
- 3、如文档侵犯您的权益，请联系客服反馈,我们会尽快为您处理(人工客服工作时间：9:00-18:30)。
A.P. Serebrov*, A.K. Fomin
Petersburg Nuclear Physics Institute, Russian Academy of Sciences, Gatchina, Leningrad District, 188300, Russia
Keywords: ultracold neutrons; neutron lifetime
2
Introduction Recently the new experiment for the measurement of neutron lifetime [1] has been carried out with a high accuracy (0.8 s). The result of this experiment 878.5 ± 0.8 s is different from the world average value 885.7 ± 0.8 s presented in PDG 2006. This difference is considerable, 7.2 s or 6.5 standard deviations. The new experiment [1] with a gravitational trap and a low temperature fomblin oil coating of the trap walls has a few advantages in comparison with the previous experiments. First of all it is a low loss factor, 210-6 per collision of UCN with trap walls. As a result the probability of losses was about 1% in comparison with the probability of neutron -decay. Therefore the measurement of neutron lifetime was almost direct, the extrapolation from the best storage time to neutron lifetime was about 5 s only. In these conditions it is practically impossible to obtain a systematical error of about 7 s. The systematical error of the experimental result [1] was 0.3 s. In connection with this situation we decided to analyze the previous experiments in more details. Most of the previous experiments with UCN storage have been carried out using fomblin oil coating at the room temperature in comparison with the experiment [1] which has been carried out at the temperature of ~120 K. The loss factor was about 30 times higher than in the experiment [1]. Very soon after these experiments [2-5] it was observed [6,7] that there is a quasi-elastic scattering of UCN in the reflection from the fomblin surface, which is rather extensive at the room temperature and can be suppressed at a low temperature [7]. This quasi-elastic scattering is happening due to surface waves of the liquid as it was proposed in the work [8]. The theoretical statements of this work have been checked experimentally in the laser experiments [8] and in the UCN experiment [7]. There is a reasonable agreement between theoretical calculations and experimental results devoted to the observation of low energy heating of UCN during the storage process in the trap with fomblin oil coating [7]. Therefore we decided to use the theory of quasi-elastic scattering process and realize Monte Carlo simulation of the first experiment [2] with the fomblin oil coating.
3
ห้องสมุดไป่ตู้
Neutron lifetime experiment [2] scheme and results Below we reproduce a short description and results of the experiment [2]. The scheme of this experiment is shown in Fig. 1. The UCN storage volume is a rectangular box, with constant height =30 cm and width =40 cm but variable length x< 55 cm. The side walls and the roof of the box are made of 5-mm float-glass plates. The oil spray head is mounted on the metal base plate and the assembly is immersed in a 1-mm-deep lake of oil. The movable rear wall, composed of two glass plates with a 1-mm oil-filled gap in between, has a 0.1-mm play with respect to the neighboring walls, except for the base plate where it dips into the oil. The surface of the rear wall was covered with 2mm-deep, 2-mm-wide sinusoidal corrugations. For half the surface the wave crests were horizontal, and for the other half vertical. This arrangement transforms within a few seconds the forwardly directed incoming neutron flux into the isotropic distribution essential for the validity of the mean-free-path formula =4V/S. The UCN inlet and outlet shutters situated 8 cm above floor level are sliding glass plates with 65-mm holes matching holes in the front wall (Fig. 1). In more details the description of experiment can be found in the work [2].
*
Corresponding author: A.P. Serebrov
A.P. Serebrov Petersburg Nuclear Physics Institute Gatchina, Leningrad district 188300 Russia Telephone: +7 81371 46001 Fax: +7 81371 30072 E-mail: serebrov@pnpi.spb.ru
1
Abstract We are discussing the present situation with neutron lifetime measurements. There is a serious discrepancy between the previous experiments and the recent precise experiment [1]. The possible reason of the discrepancy can be connected with a quasielastic scattering of UCN on the surface of liquid fomblin which was used for most of the previous experiments. The Monte Carlo simulation of one of the previous experiments [2] shows that the result of this experiment [2] has to be corrected and instead of the previous result 887.6 ± 3 s the new result 880.4 ± 3 s has to be claimed.