一辈子结束的那一天

2005年7月31日写完，2005年8月2日修改

　　自从那天在国家实验室昏倒被紧急送入医院，叶老已经在特护病房中住了十来天了。为了建设这个国家实验室，叶老呕心沥血，没日没夜的工作了近十年，国家实验室即将建成时，他却进了医院。叶老的病来得很猛，医生竭力抢救，病情却一直无法稳定。叶老时而清醒，时而昏迷。一旦醒过来，他总会问问国家实验室首次实验准备得如何了。
　　这天是叶老入院的第十四天，他的病情有所好转。在叶老的强烈要求下，医生终于允许他听取实验的准备报告。得知实验并没有因为自己的病倒而受影响，准备工作一切正常后，叶老很高兴，精气神也上来了。他又通过可视电话与实验的几位主要技术负责人通话，交待几项他临时想到的可能出现的疏漏，要求他们确保万无一失。叶老始终是病人，身子虚弱，这段精神头过去之后，他又昏睡过去。恍惚中，叶老感到自己似乎回到了六十年前，那时叶老还在上初二，别人也不叫他叶老，而是叫他……
　　“小东，小东……”，正趴在课桌上睡觉的叶小东被同桌吵醒，正准备K他一顿，突然发现自己被物理老师点起来回答问题。望着这个刚刚大学毕业分配来的年轻老师，叶小东决定赌一把。
　　“老师，刚才那个问题我没有听清楚，您能不能再说一遍？”
　　“行，我刚才问的是请举一个例子，说明一下由于参照系选取的不同，运动也不一样。”叶小东赌赢了，这位新老师没有象语文英语老师那样奚落学生一番，而是按照叶小东的要求复述了一遍问题。
　　刚刚从睡梦中醒来的叶小东听完问题后，仍旧感到昏昏呼呼的，于是猛的甩了甩头，又是一阵天旋地转。看着旋转的天地，叶小东有了主意。
　　“老师，比如我在这里旋转，”说着，叶小东做旋转状，“旁人看起来，我在转圈，”说得兴起，叶小东真的转了起来，“可是在我眼里，是整个世界在旋转。”这句话刚刚说完，一不小心拌了一下，给摔到地上。从叶小东开始转圈就忍着不敢笑的同学们终于忍不住都笑了起来。
　　叶小东感到很没面子，只好拍拍身上的灰尘，灰溜溜的站了起来，准备接受老师的批评。他自己也觉得刚刚那个回答太异想天开了。出乎所有人的意料之外，物理老师先是询问了一下叶小东有没有受伤，然后把他大大的表扬了一番。叶小东不禁有点飘飘然，难道刚才的那一跤摔得老师大发慈悲不再计较自己上课睡觉没有答对问题？没有人注意到，当叶小东说出答案后，物理老师眼中那惊叹的目光。
　　既然老师这么给自己面子，当众表扬自己，叶小东决定也给他面子，剩下这半节课认真听他上课。叶小东今年12岁，刚上初二，这是他有生以来第一次接触物理学。叶小东感到黑板上的那些物理定律似乎有魔力似的，把自己牢牢吸引住了。半节课很快就过去了，刺耳的下课铃声把叶小东从物理世界中拉了回来。回过神之后，他突然想到物理课是下午的最后一门课，赶快闪人，不然就没位子了。叶小东把文具盒书本往书包中一塞，一把抄起书包，斜背在肩上，冲出后门，直奔校外的游戏厅。
　　20世纪九十年代初，在中国无数的城市中，一夜间冒出了许多游戏厅，“街霸”，“三国志”等游戏风靡万千少年，叶小东就是其中一员。远远望见游戏厅门前攒动的人头，叶小东就知道已经晚了。不过他还是走入游戏厅，希望能够蹭到一台机子过过瘾。游戏厅内，“Hogan”之声不绝于耳。一听到这个，叶小东就从心里痒痒。转了一圈，找不到一台空机子，看样子一时半刻也不会有人走。叶小东决定转移战场，希望另外那家游戏厅能够有空的机子。
　　刚刚走出游戏厅，就发现物理老师往这里走来，叶小东低低的喊了一声：“王老师好！”然后站在老师面前低下头，准备挨训。
　　“你跟我来。”没有预料中的训斥，只是这么一句，不咸不淡，然后王
老师转身向学校走去。叶小东忐忑不安的跟在后面，心中直打鼓。很显然，王老师的目的地是物理教研室。走进教研室，其他老师都下班回家了，偌大的房子里，只有王老师和叶小东两个人。
　　王老师径直走到自己的办公桌旁，翻看着什么东西，拿着教案回过头来，发现叶小东在自己面前低着头站着，一副丢了魂的样子。他装作没有看到，指着一把空椅子说：“小东，怎么还站着，坐吧。”叶小东坐下后，王老师接着问：“怎么样，今天上课听懂了么？我现在来考考你。”然后就考起叶小东物理知识来了。
　　叶小东心说，这老师真是奇怪，怎么抓到我玩街机，就像什么事也没有发生似的，幸好今天上课后半节课我听了，不然就不好过关了。然后老师问什么，就照课堂上留下的印象回答。突然，王老师问到：“课堂上我问的那个问题，你是怎么想到那个答案的？”叶小东见周围没有旁人，就一五一十的把当时的情形告诉了王老师。
　　听完叶小东的解释，王老师哈哈大笑：“小东，你在物理上很有天赋，大好时光，用来玩电子游戏太浪费了。你以后每周二，四物理课下课后不要回家，到我宿舍来，我来教你更好玩的东西。”然后摆摆手，放叶小东回家了。
　　走出办公室，叶小东长吁一口气，游戏厅是不敢再去了，还是赶紧回家吧。叶小东的家在市棉纺厂。由于美国的制裁，上个世纪九十年代初，全国棉纺织企业都亏损，叶小东父母所在的棉纺厂也不例外。厂里发不出工资，作为技术骨干的父亲只好出去四处找活干，母亲守在家里，却越来越管不住渐渐长大的小东。聪明伶俐的叶小东小学时成绩优秀，尤其是数学和自然课非常好。为了孩子学业，父母找了不少关系，花了不少钱把叶小东送到了现在这所重点中学念书。
　　刚入初中时，一直成绩优秀的叶小东在摸底考试中考得很差劲。所谓强中更有强中手，在全市成绩最优秀的学生中间，叶小东只属于末游。让父母失望的叶小东自己也有点自暴自弃。在班上几个哥们的鼓动下，叶小东学会了玩街机。被生活所迫压得喘不过气来的父母顾不了小东，他也就一发不可收拾的沉迷于街机中。
　　叶小东回到家，却没有见到母亲。母亲在桌子上留了张纸条，要小东自己自己先吃饭。唉，看来又被人叫去打麻将了。叶小东只好从电饭煲中拿出热好的剩饭剩菜，打开电视机，边吃饭边等着六点半的动画片《天空战记》。
　　吃完饭，动画片也看完了，天色渐晚。“真无聊！”叶小东嘀咕着。突然想到后天星期四物理老师要找他，不禁又犯慌。算了，看书去吧，叶小东心中说到。
　　晚上十点，叶母赢了一把“清一色”后突然记起孩子还在家没人照顾，于是不顾牌友的再三阻拦，离席回家了。叶母进门后发现孩子房中灯还亮着，推开房门，走到小东的身后才发现小东正在温习物理，心说：唉呀，这孩子终于开窍了，懂得自觉读书了。叶小东这才发现母亲已经回来了，正要起身，叶母赶紧按住他说： “别动，好东东，继续看书，我给你倒水去！”叶小东看看钟，才发现已经看了三个小时的书了。
　　那天晚上，叶母不断的夸小东开窍了，搞得叶小东怪不好意思的。小东晚上睡觉时，做了一个甜甜美美的好梦。
　　叶老从沉睡中醒来，眼前是特护病房的输液瓶和洁白的墙壁。朝阳升起，阳光射入病房，在灿烂的阳光照射下，新的一天开始了。
　　叶老突然心有所动，问了问正在整理病房的护士：“今天几号？”当得知已经进入7月份时，他知道距离预定7月10日开始的实验室首次实验没有几天了。叶老还能够清楚的记得，59年前的七月，当时的叶小东，第一次体会到了物理实验探索的魅力。
　　自从叶小东被王老师要求每周二、四去开小灶补习后，他的物理成绩突飞猛进，顺带的数学，英语，语文等诸科成绩都大有进步。这次初二下学期的期末考试中，叶小东更是进入了全班前十五名。叶母从儿子那里得知有这么一位认真负责的物理老师后，非常感激，硬是拉着儿子找到物理老师，准备请老师吃顿饭以作答谢。王老师抵不住叶母的盛情，只得吃了这顿饭。
　　在饭桌上，王老师对叶母说：“你儿子有天赋，应该在物理方面继续发展，既然我吃了这顿饭，我就要负这个责。暑假你把儿子交给我，我专门来教他。”叶母一听，有这么好的事，自然是答应了。
　　再说叶小东，自从跟着王老师学习物理后，游戏厅是很少去了，原因在于做起老师留的物理习题就经常忘了时间。这次期末考试成绩这么好，他正和几个哥们约好了准备暑假好好去游戏厅过过瘾呢。一看王老师又来了这么一个提议，虽有一千个不愿意，当着老师和母亲的面，只得应承下来。于是叶小东每周一到周五有五天白天都得到王老师家去。
　　第二天正好是周一，叶小东一大早就被母亲叫醒，让赶快去王老师家。敲开王老师的家门，王老师见是叶小东，说：“小东来了，你先做暑假作业。”然后打开箱子捣鼓着什么。叶小东很好奇，想偷偷看两眼老师在干什么，却被老师的身体给挡住了。只好从书包中拿出物理暑假作业，心不在焉的做着。
　　一个钟头过去了，王老师终于停了下来。“来，小东，你过来看一下。”叶小东巴不得看看王老师在干什么。咦，怎么是个风车一样的东西，不同的是叶片是金属的。王老师转起了“风车”，很快，风车片开始噼里啪啦的冒起火花。“小东，站在这双胶鞋上，用手摸一下这个金属球。”叶小东看到 “风车”连在一起的金属球，畏畏缩缩的不敢走上前来。
　　“怎么，怕了？那好，你来摇动静电发生器，我来试试。”王老师站起身来，示意叶小东过来继续摇。叶小东摇动了静电发生器，看着那时不时闪现的火花，他感到既刺激又好玩。王老师踩在胶鞋上，手指轻轻的一触金属球，只见王老师那近两个月没有理的头发都竖了起来。“哇！”叶小东惊叹起来。
　　“王老师，你的头发怎么都竖起来了？”叶小东问。
　　“因为静电。我手触摸金属球头发就带上了静电，带相同电荷的头发丝之间相互排斥，就竖起来了。”接着，王老师详细的向叶小东讲解电的性质。听着王老师的讲解，叶小东感到脑子都有点不好使了。他时不时地问几个幼稚的问题，王老师也不笑话，而是耐心的解答问题。
　　那个暑假，真是一段梦幻般的日子。王老师带着叶小东遨游在物理科学的奇妙世界中。王老师简直是个魔术师，而他那十来平方米的简陋的宿舍就是一个魔术屋，不论要做什么实验，似乎都能从堆在宿舍的几个大箱子中找到合适的工具实现它。此外王老师还告诉他，许多人们习以为常的事情中都有不寻常的物理规律在起作用。天为什么是蓝色的，星星为什么会闪烁，为什么会有打雷闪电，寻找这些为什么，回答这些为什么，就是物理学要做得事情。
　　叶老还记得有一天晚上，王老师没有让他早早的回家，而是带着他去了街心公园，在草地上躺了下来。王老师指着天上的星星说：“看，天上这些闪闪发光的星星，叫恒星，它们都是远方的太阳。不论地面上人事如何变幻，星空是永恒的。我相信这外部世界变化的规律也是永恒的，令人心醉的，值得用一生来探究。”叶小东似懂非懂的点了点头，顺着王老师的指尖望去，似乎有某种神圣的东西从星空中通过双眼传入自己的身体中。
　　暑假后，叶小东彻底的迷上了物理。从王老师那里了解到要学好物理，数学，语文，英语都得学好。因为物理规律是用数学语言描述的，而最新的物理，大都是用英语书写的。于是除了物理，叶小东也努力钻研其它各门课。不得不说王老师独具慧眼，在叶小东的刻苦学习之下，潜力都发挥出来了。叶小东的成绩步步攀升，很快稳定在班级前五名。
　　一个偶然的机会，叶小东从老师谈话中了解到，王老师的来历大不简单。王老师是近十年前这所初中毕业的高材生，高中毕业后考入了首都的那所理科名校学物理，92年大学毕业时不知犯了什么错误给打回原籍。他在家待业一年准备考研。考研失败后，最终还是母校收留了他。叶小东猜测了好久始终猜不到王老师犯了什么错误。
　　回忆到这里叶老心中一阵抽搐，护士也注意到仪器显示的指标出现异常。
　　“叶老，您怎么了？”
　　“没事儿，没事儿。”叶老摆摆手，从回忆中醒来。
　　“叶老，您病还没好，不要总记挂着工作。安心养病最重要！”护士见叶老很疲惫，连忙开导他。
　　叶老沉默了一会，又说：“小张啊，把主治大夫叫来，我有事儿想要问他。”
　　张护士连忙找来正在值班的主治大夫。叶老双眼直直地望着主治大夫，半晌才开口：“医生，我到底得了什么病？”
　　医生连忙说：“叶老，您别胡思乱想，安心养病。”
　　“我明白，你们有规定不得向病人透露病情。可是我心里明白是怎么回事。人老了，岁月不饶人，我身体的零件快报废了。哪怕现在医学已经克服了癌症和艾滋病，可是它对抗不了衰老。”
　　医生欲言又止，叶老接着说：“我知道自己大限将至，我只希望医生你能帮助我活过7月10号国家实验室首次运行的那一天，见到自己为之奋斗了一辈子的东西变成现实。如果能那样，我将了无遗憾。”
　　说完这些话，叶老的精力似乎用尽了，头一歪，又陷入了昏迷。医生见状，马上叫来护士抢救。
　　叶老正在步入他人生最后的一段旅程，58年前的叶小东，在他生命中头一遭，亲眼目睹了“死亡”。
　　自从叶小东得知王老师的身世后，就特别留意他的一举一动，很快发现王老师的特意独行。王老师很少跟其它老师打招呼，下班成天猫在的宿舍不出来走动，而且除了物理和教学他几乎不关心任何其它事情。叶小东暗地里也很奇怪，王老师人长得帅气，个子高，又那么有才，可怎么工作了好几年没见王老师找个女朋友呢？一次偶然的机会，叶小东去王老师家时发现王老师枕头旁有一张照片，照片上的女孩清秀脱俗，背景是天安门毛主席像。叶小东明白了，原来王老师是有恋人的。
　　初三的复习是很紧张的，叶小东很少有机会再去王老师家做物理实验了。不过他已经计划好了，等中考一结束，就去找王老师，试着把去年做了一半的半导体收音机给做出来。
　　中考结束了，成绩出来了，叶小东顺利的考上了省重点高中。他兴高采烈的去找王老师报喜，却发现王老师不在。向其他老师打听后才知道王老师去北京了。叶小东估摸着他肯定是会女朋友去了。
　　没有找到王老师，叶小东有点失落。不过这两年叶父在外边打工赚了不少钱，他家的情况好了许多。父母见小东考得这么好，于是决定请几天假，带着小东去北京旅游。在北京，叶小东特意来到了王老师曾读过书的那所高校，在那里他暗暗发誓，以后也要到这所学校读书。
　　叶小东本来期望能够在北京巧遇王老师的，可是北京这么大，直到离开北京回家还是没有遇到王老师。回到家后，叶父要小东赶快去王老师家看看王老师回来没有，回来的话跟王老师说下周一庆祝小东考上重点高中请客，要王老师一定到。
　　叶小东不会忘记，那一天是1995年7月10日。他是吃了晚饭离家去王老师家的。在王老师门口敲了半天，没有人答应，正准备走时，门内传来王老师的声音：“谁呀！”，叶小东连忙说：“是我呀，小东。”叶老师打开了门，把叶小东让到屋内。
　　叶小东坐下后，把父亲的话复述了一遍，王老师只是嗯了一声，却没有答应。这时屋内的气氛有点诡异。叶小东仔细看了看王老师，从北京回来后，王老师整个人都瘦了一圈，眼睛周围是大大的黑眼圈，似乎好几天没有睡觉了。叶小东猜测王老师可能是失恋了。见王老师心情不好，叶小东准备马上开溜。
　　“王老师，就这么说定了，到时候我来接您。您刚从北京回来，肯定累了，早点休息吧。”说完，叶小东就往外走，谁知却被王老师给叫住了。
　　“小东，你别走，陪我说说话。”王老师的话语中透露着一股辛酸。叶小东停住了脚步，走回来坐了下来。
　　“小东，你这次考上重点高中，值得庆祝。可是不能松劲，要继续努力呀。当年我也是从这所初中考上你即将去的那所重点高中，然后考上大学的。”说到这里，王老师顿了一下，“我猜你已经听说了，我考上了名牌大学物理系，毕业后却被打回原籍，回到这里教书。我不甘心呐！”说到这里，王老师已经带着哭腔了，然后强行忍住，也许是感到在学生面前失态了，王老师接着说：“小东，你太单纯了，象我当年一样，这样会吃亏的。你以后对人对事要多留个心眼。不管怎么说，能够在这里碰到你这样一位真心热爱物理的学生，是我的幸福。我的物理之路走到了尽头，可你的还刚刚开始。我猜你已经下决心以后报考物理专业，对么？”
　　叶小东点点头。王老师欣慰的说：“看来我没白给你开小灶。你回去跟你父亲说，下周一我就不去了，我教你是老师的责任。”叶小东还想再劝一劝王老师，却被王老师二话不说，送出了家门。
　　从王老师家中走出来，体会这老师的那一番话，叶小东总感觉老师那温暖的话语中有些怪怪的味道。刚刚出楼道，叶小东突然感到肚子痛，估计晚饭吃坏肚子了，情急之下从学校阅报栏中撕下几张通知就冲入20米外的公共厕所。
　　半个小时后叶小东才从厕所中出来。夏夜的校园，教学楼一片漆黑，路灯昏黄，只照亮了底下几米见方的区域。没来由的，叶小东感到一股寒意，这么大的校园就我一个人，唔！还是赶快回家吧。
　　这时，远处王老师家大开的窗户引起了叶小东的主意。王老师家灯全亮了，突然，王老师出现在窗台上，几乎没有犹豫，从五楼的窗口跳了下来。叶小东惊呆了，眼见着王老师落到二楼时，突然像是有一个洞凭空出现，王老师消失在空气中，只有一阵狂风夹着闪电继续吹到地面，将楼底下的尘土高高扬起，落叶也被点燃，吹到空中。
　　叶小东呆呆的站了几分钟才明白过来，王老师跳楼了，却又莫名其妙的消失了。后来这段日子是混乱的。叶小东报警后，没有人相信他的“幻觉”，但是公安干警确实在王老师的宿舍找到了一封遗书，上面说自己走后，留下的书籍和实验材料送给叶小东。警察分析，王老师是在得知恋人离开自己傍大款之后收不了打击，所以才会自杀的。可是警察始终却找不倒王老师的尸体，最后只好按照失踪案处理。依照王老师的遗书，王老师的家属把他留下的物理书和物理实验器材送给了叶小东。
　　尽管警方得出了王老师的自杀原因，可是叶小东后来自己分析，这只是原因之一。在警察调查王老师失踪案时，叶小东才知道为什么毕业于名牌大学的王老师会被打回原籍。临近毕业时，成绩优秀，热爱物理的王老师本来已经被系里选中保送读研的，可是却被人告密说他曾经参与过那场政治风波，政治上有问题。学校不理王老师辩解说那天只是好奇去看了看，直接将他的档案打回原籍，并把这件事记录在档案中。这也是后来王老师考研落榜的原因。叶小东认为，失去在物理上继续学习的机会后，又丢掉了爱情，双重打击令王老师失去了对人生的希望。叶小东甚至想，王老师见到自己能够走上物理这条路，是他死前最后一点慰藉。
　　多年以后，叶小东意识到，自己的一生就是在1995年7月10日结束的。那一天，他在王老师面前立志一辈子从事物理，他找到了自己生存的目的，于是他的一生也就结束了。
　　叶小东在物理学上的发展并非一帆风顺。由于受到王老师自杀的影响，叶小东没有考上王老师的母校，而是来到了西北的一所高校学物理。在大学中，叶小东一次次的回忆当年自己见到的那一幕，坚信那不是幻觉。后来他偶尔从互联网上发现了爱因斯坦和Rosen曾经构造过“虫洞”这么一个模型，似乎能够解释王老师为什么消失。叶小东疯狂的收集这方面的资料，并决心以后从事“虫洞”研究。
　　本科毕业后叶小东考上了中科院理论物理所的直博生，研究广义相对论。上博士以后，面对课题中的许多困难，叶小东又一次回想起王老师第一次夸他在物理方面有过人的天赋的情形，他不禁一笑：那时的王老师不可能仅仅通过一个问题就能断定自己是天才，那不过是作为一个老师给学生的鼓励。后来自己被王老师激发出学习物理的热情后，王老师才会主动提出暑假时继续给自己补习物理，那时王老师才真的认为自己适合学物理吧。想清这一点后，叶小东没有气馁，而是更加投入到研究中去。叶小东博士期间做得相当好，在广义相对论和宇宙论方面作了不少出色的工作。
　　博士毕业叶小东他去了意大利国际理论物理中心做量子引力方面的研究，做完两年博士后，叶小东回到母校做副教授，继续研究量子引力。21世纪初，物理学与各个学科交叉非常多，做量子物理的许多都转去做量子信息，量子计算，量子通信，而叶小东一直守着量子引力这一块儿，研究“虫洞”。不过这个领域似乎已经被发掘完了，回国后叶小东一直没有发表出色的论文。
　　由于回国后论文发得很少，叶小东第一次评教授失败了。那时他也老大不小了，却一直没有结婚，他的心全被“虫洞”给迷住了。后来扭不过父母，回家相亲后，跟家乡的一个女孩结婚了，那时叶小东已经35岁了。妻子很贤惠，很宽容，从不干涉叶小东的工作。
　　新婚后的叶小东似乎有了动力，多年研究虫洞的积累在这时解出了硕果。他写出了一篇突破性的论文，指出利用现有的实验条件，人们已经可以制造控制虫洞了。他的论文证明了有另外一种虫洞模型存在，能够有效的克服原有模型中虫洞内需要负能量才能稳定的问题，他把这种虫洞称之为“双向虫洞”。在这个“双向虫洞”模型中，虫洞联接的两个时空点间的能量质量交换是双向的，且两点各自的局域能量与动量密度在交换时守恒。在后来的一篇论文中他又发现，人工制造出来的虫洞，另外一端的时空坐标是不确定的：也许在宇宙的另外一头，也许相隔万年。可惜当时物理学界并没有多少人注意叶小东的这几篇论文。不过能够在物理学顶级刊物上连续发表两篇论文，已经足够叶小东成功的评上教授了，这时他已经四十多了。
　　日子就这么平平淡淡的过着的：叶小东继续教他的书，偶尔发表几篇关于量子引力的论文，参加一下学术会议，慢慢的在学术界也小有名气。于是他开始尝试申请基金支持建设实验室来验证自己的理论。可是国内没有人相信按照叶小东的设计，真的可以造出虫洞来。很多人劝他：转个方向吧，小东，物理学这么广阔，能够研究的东西多着呢。不得已叶小东转到了一个相近的领域，可是他从未放弃研究虫洞，一直在暗地里跟踪相关进展。十几年后，叶小东近六十岁时，美国的一个科研小组宣布，依照叶小东当年的论文，他们成功的制造出了一个直径1mm的双向虫洞，叶小东的预言完全实现了。
　　第二年的诺贝尔物理学奖授予了叶小东和那位首先制出双向虫洞的美国物理学家。于是叶小东成为了21世纪首位获得诺贝尔物理学奖的中国人。叶小东年纪大了，地位又高，慢慢的对他的称呼也变了，年轻一辈的开始叫他“叶老”了。他知道，当被人尊为“老”时，自己有创造力的学术生命也就到尽头了。
　　然后，叶老开始利用获得诺贝尔奖的机会申请到了足够的资助开始建实验室，自己造“虫洞”。叶老知道，这是自己一生中最重要的一个机会，也是最后的一个机会，让他再见一次“虫洞”。这个秘密他藏在心里，谁也没有告诉。
　　知易行难，尽管叶老创建了理论，可是他毕竟不是实验专家。为了建设实验室，他不得不学习做一个实验物理学家，此外还得学习如何调度安排人力，最高效的建设实验室。不同于美国科学家的那个实验，他们现在造的实验装置预定要产生直径一米的“虫洞”，实验需要的能量要高出若干个数量级。叶老他们只有先重复造出1mm的小虫洞，然后一方面提高能量，一方面提高能量利用效率。积累的近十年，才制造出现在这个实验装置。
　　叶老昏迷了3天，在医生的抢救下终于恢复了意识，但这次他的精神比上次恢复意识时差多了。
　　“几号了？”叶老用嘶哑的嗓子问护士，也许是他的声音太低，护士把头低下来，凑到叶老的嘴边，要求叶老再说一遍。这次，护士终于听清了。
　　“10号了，叶老。”护士回到说。
　　“啊，已经十号了，实验是不是已经开始了？”叶老声音中有点焦急。
　　“还没有，离预定时间还有一个钟头。”听到这个消息，叶老放松了下来。然后要求护士接入实验现场的摄像信号，从全息电视上观看实验。
　　国家实验室现场，已经进入倒计时阶段。负责指挥调度的实验室副主任是叶老的学生，十多年来一直跟随叶老。在叶老入院后他心急如焚，却由于工作需要脱不开身前去探望。
　　时间一分一秒的过去，实验装置开始启动，预热，输入初始条件，在计算机的控制下，具有太电子伏特量能量的粒子绕着环路奔跑，越跑越快，最后在终点处，两束粒子相碰，碰撞处发出强光，然后一切归于平静。可是探测器却把这一瞬的情形记录了下来。在碰撞的一瞬间，空间打开了一个直径约一米的虫洞，持续时间大约1毫秒。由于“叶氏不确定性定理”，另外一端的时空坐标可能取值将遍及整个宇宙，直到虫洞形成的那一刻才会塌缩到某特定的点上。在这个过程中，换算成国际单位制，有大约70kg的质量能量发生了交换。
　　实验结果是在实验结束后两个钟头由计算机处理得出的。实验室副主任已经知道叶老正在病房观看实验，于是决定把结果通过大屏幕显现出来。
　　正在观看实验的叶老看到了大屏幕上的实验结果。猛然间他似乎想到了什么。70kg，按照王老师的身高体型，大约就是这么重。难道自己当时看到的虫洞就是现在这个？通过虫洞的极端高温高压环境后，王老师在一瞬间分解成基本粒子，而那一阵风与闪电，很像是通过虫洞由实验室交换到那里的高能粒子流。
　　想到这里，叶老几乎可以肯定了：他花了10年工夫建设的国家实验室成为了他物理启蒙老师的墓地。可惜他已经没有时间验证他的推论了。于是他笑了起来，笑他自己，笑这可恶的命运。他又一次的回顾自己的这一生，回忆起那个梦幻般的暑假，回忆起与王老师相处的那段日子。他不知道，是自己杀死了王老师，还是王老师自杀。原来， 1995年7月10日，真的是我的最后一天，这是叶小东最后的意识。他带着轻蔑的微笑，陷入了永恒的沉眠。
(全文完)

我的第一个研究工作发表了

2005年上半年，我在导师指导下完成了第一个工作。我们本来打算在那个暑假后整理投稿，但在这个工作的启发下，我很快又发现了更加值得做的题目，于是才有了我的第一篇PRA。再后来就是我申请出国，这个工作也就彻底放下了。有关这段事情，我以前也记录过。

到了2009年，我导师让师弟整理了一番，投稿，最近终于发表了。从接收到发表差不多2年。而从我做这个工作，到现在已经6年多了。如果读了这篇论文，和我的第一篇PRA比较，就会发现二者是有很密切的联系的。我之所以能够想到用多原子系统完成量子计算可以起到加速作用，就是受了这个工作的启发。

时空隐身斗篷

（首发于“果壳网”）

当风靡全球的电影“哈利•波特”中魔法师的隐身斗篷头一次出现在电影荧幕时，人们还以为这是魔法和幻想。可如今随着超颖材料（metamateria）技术的发展，这已经不再是幻想了。人们已经可以在微波波段设计出隐身斗篷，使得电磁波从障碍物旁绕过去，达到隐身的效果。随着纳米技术的发展，制造出能够在光波波段工作的三维隐身斗篷也不再是幻想了。对电磁波隐身的斗篷工作原理如图所示。

/gkimage/du/t9/va/dut9va.png

如图，光波沿着平行于x轴的方向传播，在0点处有一个隐身斗篷。当光传播到0点附近时，其传播路径发生的弯曲，贴着斗篷的表面传播。当光远离斗篷时，又恢复为平面波。这样在远处的观察者就根本无法发觉原点处的隐身斗篷以及其中我们要隐藏的物体。

我们知道时间与空间是相关联的。既然我们可以制造出隐身斗篷把空间上某点隐身起来，那么是不是也可以制造出一个时空隐身斗篷，用于隐藏发生在特定时空点上的事件呢[1]？如果可以那真是太酷了！想象一下，这是任何推理小说和特工电影都没有呈现过的景象：在一个受到严密监控的密室中，找到一种办法把某个将要发生的事件屏蔽起来，使它不被监视器记录。这就相当于偷偷地进入那个密室，打开监视器，清除相关内容，然后关门离开，这样在监视器中就像是什么也没有发生一样。只不过一个是事前屏蔽，一个是事后消除。

这听起来就像是科幻小说一样，可实际上用同样的超颖材料技术，从原理上来说，不是不可能实现。科学家们把它称之为“时空隐身斗篷”，或者说是“历史编辑器”。最近的一个实验也从原理上验证了时空隐身斗篷的可行性[2]。其基本的原理其实很简单。与传统的空间隐身斗篷不同，时空隐身斗篷并不弯曲光纤在空间中的轨迹，它利用的是超颖材料来控制光，制造出暂时的时间空白点来隐藏某个发生了的事件。具体来说，我们可以把探测的光分为波前领头部分，它被加速了，于是在事件发生前就跑过去了；还有波尾部分，它被减速了，在事件发生后才通过那个区域。最终把波前和波尾合并起来，就抹去了时空隐身斗篷存在的痕迹。在一个远处的观察者看来，任何在时空斗篷存在时间内存在的事物信息都被抹掉了。换句话说，对那个特定的时空点，我们创造了一个隐身斗篷，其中任何能量、信息和物质的痕迹都无法被观察到。我们也可以用下图来描述这个过程。

/gkimage/50/0d/9g/500d9g.png

上图，光在时空中传播，横坐标是x，纵坐标是时间轴ct，c为真空中的光速。当光达到x=0和ct=0附近时，就发生了时空上的弯曲，绕过了原点。换句话说，发生在时空原点附近的事件被这个时空斗篷隐藏起来了，扭曲了我们所观测到的物理现实。

目前的实验只能实现一个空间维度加上一个时间维度的时空隐身斗篷，未来也许可以实现三维空间+一维事件的隐身斗篷，这说不定将打开科幻的新领域。设想一下，有一天这种时空隐身斗篷隐藏在我们周围，那么依赖眼睛来观察世界的我们，所看到的将不再是原来世界，而是被扭曲和编辑过的世界。又或者说宇宙中本身就存在类似的这种时空隐身斗篷，扭曲了天文学家所看到的天幕，让我们看到的都是假象（像不像“三体”系列里的智子干的事呢）。目前的时空隐身斗篷只抹去了光所携带的事件信息，所以还可以用其他媒介手段来消除隐身斗篷的影响。可我们也不能否认存在其他类型的隐身斗篷，把引力波，中微子等诸多媒介所携带的信息都抹掉。实际上，我们自己已经做出了声波的隐身斗篷，未来有一天也许可以用于抵御海啸[3]。

参考文献:

[1]	A spacetime cloak, or a history editor，Martin W McCall et al, 2011 J. Opt. 13 024003， http://iopscience.iop.org/2040-8986/13/2/024003/

[2]	Demonstration of temporal cloaking，Moti Fridman et al. http://arxiv.org/abs/1107.2062

[3]	Broadband Acoustic Cloak for Ultrasound Waves, Shu Zhang, Chunguang Xia, Nicholas Fang, Phys. Rev. Lett. 106, 024301 (2011)， http://arxiv.org/abs/1009.3310

读博士的体会

在本文中我将与大家分享一下读博士的经验和体会。包括几个部分：首先，为什么读博士；其次，读博士的注意事项；最后，读博士的收获。

在读博士之前，应该明确自己的动机，并确信自己知道读博士付出与回报是值得的，并愿意承受种种付出的代价。作为理工科的学生，我读博士的动机是希望能够更加深入的学习和了解自己的爱好：物理。而且也只有读完博士，经过了严格的训练之后，我才能在相关的研究机构里面找到合适的职位，从事物理学研究。我也知道，念博士的代价是让自己更晚的离开校园，收入很少，而且也缺乏结交异性的机会。但是为了实现自己的理想，我愿意承受这些代价。在你念博士之前，要想清楚自己要什么，为了实现自己的理想，又要付出什么代价。权衡之后，再做出决定。

读博士的注意事项，或者说一些技巧，这部分是重点内容。所谓工欲善其事，必先利其器，在念博士前了解一下相关的小窍门，会起到事半功倍的效果。首先，我们要明确，博士教育的培养目标是一个具有独立研究能力的学者，所有的教育和培训项目都是围绕这个目标而来的。而所谓的教授，就是已经用自己的独立研究能力做出了很多成果的人。可见，博士与教授也就是一步之遥，欠缺的也就是火候。在念博士时，你是你导师的研究助手和合作伙伴，导师布置的任务是很好的研究训练，一定要主动的、及时的、认真的完成。很多人的弱点是，主动性不够，害怕在导师面前暴露自己的不足。在念研究生时，与导师的沟通不够。念博士时，必须要特意的来训练自己，主动的积极的找导师，与导师讨论。据我观察，导师都是喜欢主动学习的学生的。不管你的水平在导师眼中是好是差，他都不会在乎，他在乎的是你是否努力，是否在不断的进步，为他的课题做出贡献。

第二个非常重要的窍门在于，你要跟谁念博士，就决定了你的学术起点。导师有很多类型的：有的导师不管学生，任学生自由发展，但同时也会给学生提供足够的研究条件；有的导师对学生非常负责任，经常督促学生，检查研究进度；有的导师学术水平一般，却又很喜欢批评学生，对学生的工作不懂，却又装懂，等等。从学术水平来说，从院士级别到普通教授相差极大。一般而言，跟着名教授，或者院士，你的学术研究起点将会很高，也会很容易进入学术圈中，给自己未来的发展提供了非常大的便利。比如说，跟随名教授念博士，你将有机会参加学界的高水平会议，甚至在会议上做报告，与业内的专家学者交流，建立自己的学术声誉。让圈内的人认识你，乃至赞赏你非常重要。因为博士毕业后，你将继续在这个圈子里面生存，利用博士期间的学术交流活动建立自己的人脉和关系，将极大的促进你的事业的发展。但同时，名教授又非常忙，他们经常没有时间指导你，所以你只能跟着师兄师姐，或者所谓的“二老板”来学习做研究。这很考验你的独立学习和研究能力。既然跟谁念博士如此重要，那么在报考博士之前，你就得很仔细的斟酌，考察各位教授的为人品性，学术水平，业内口碑等等。千万不能考那些既没有研究基金，也缺乏学术成绩，而且对学生不好的导师。

第三个窍门，是如何看待训练独立的研究能力。独立研究的能力，包括分析问题的能力，掌握信息的能力，学习新知识的能力。这些能力中，掌握信息的能力，其实就是如何检索文献，找出与你的研究相关的前人研究工作的能力。现在的网络数据库很多，大家的检索办法都大同小异，建议大家好好学习Web of Science数据库的检索办法，就可以触类旁通。我个人体会，迅速的检索出正确以及完全的文献，是完成研究课题的助推器。拿到文献之后，先要粗读，看看标题和摘要，找出感兴趣的。然后把比较感兴趣的文献的前言和结论仔细读读。注意，此时要做好文献笔记，用于以后查阅。念博士期间都是在学中做研究，以及在研究中学习的。很多高深的知识和理论，在真正做研究时可能只需要用到其中一小点，你只要把那一点给学透了，就行了。分析和解决问题的能力，得靠观察学习导师如何分析问题，师兄同学如何分析问题来慢慢体会。

第四个窍门，是积极的整理笔记，并按照论文的格式来写笔记。念博士免不了写论文，而人并不是天生会写论文的。论文的格式与要求与高中语文所教的作文的要求有很大的不同。学术论文是有约定的体例和严格的语言要求的，甚至你可以说它有八股文的味道。要想写好论文，把自己的思想严谨而规范的表达出来，必须经常训练。那么平时的研究笔记，课程学习要求的大作业等等都是训练写论文非常好的途径。每隔一段时间，把自己读过的文献笔记整理起来，按照文章的格式和要求来写，这对自己是很大的训练。我个人的体会是，做研究和课题时，写论文是最痛苦的事情，同时也是最磨砺自己思想的时候。因为很多东西，你不写下来，是不知道自己到底弄明白没有的。只有你能够把它清晰的表达出来时，这个知识才算是真正掌握了。而很多体会和研究结果，不写下来，也不知道到底研究透彻没有。在论文写作过程中，会令自己从各个角度来分析问题，思考同一个问题，直到对问题烂熟于胸。

现在，我已经博士毕业了，我觉得收获很大。念博士，首先是一种逻辑思维的训练，以及分析问题和解决问题能力的训练。我虽然学的是物理学，但我相信自己具备了独立学习和掌握其他学科研究的能力。未来如果工作需要或者自己感兴趣，我能很快的进入新的学科学习和研究。其次，我结识了一批学术同行，与志同道合的一些朋友建立了学术联系，有了自己的圈子。最后，我靠这个文凭找到了一份适合自己发展的工作，有了稳定的收入。我个人觉得，念博士后最重要的收获就是一份自信——自己能够在学术和相关领域取得成绩的自信。

后记：这篇文章是我博士毕业后不久写的，现在看来似乎还有些价值，就发到blog上跟大家分享。

我申请到基金了！

前天刚刚收到自然基金委员会的email通知，我申请的青年基金项目获得了批准，同意资助我24万元，项目时间是2012年1月到2014年12月。我能够申请到这个项目，离不开我的导师们多年来的培养和支持，也离不开我的同事们的帮助和合作。在此一并感谢大家！我会努力工作，做出好的成绩，不浪费这些纳税人的血汗钱。

Quantum Information processes in atomic ensemble

In this note I will summarize the works by Polzik’s group. I will give the detials of the theoretical treatment in their works.

1. Introduction

Atoms are good at storing and processing quantum information because of its long coherent time. Light is suitable for long-distance quantum state transform. Hence, the coupling between atomic ensemble and light is the key factor for realizing quantum network.

As the number of ensemble is huge, we can treat it continuously, just as we treat light. Because of the collective coupling, the quantum information processes, such as teleportation, can be archieved without cavity, which is greatly decreased the difficulty of the experiments.

1.1. continuous discription of light and matter

Let’s firstly discuss the light/matter interaction. We will use continuous discription of the light and matter. We quantized electromagnetic field in one dimension and only care about on polarization mode (x-polerization)

$\displaystyle E= \sum_{\lambda} \sqrt{\frac{\hbar \omega_\lambda}{2 \epsilon_0 A L}} \big( a e^{ikz} + a^\dagger e^{-ikz} \big), \ \ \ \ \ (1)$

where {}{}

A, cross seciton of the field mode;

L, length of the mode;

k, wave vector. Impose periodic boundary conditions, the ${k}$ space resolution ${\Delta k = 2\pi/L}$ . When ${L \rightarrow \infty}$ , ${\Delta k \rightarrow 0}$ . ${k}$ space tends to continuous. We make changes

$\displaystyle \sum \Delta k \rightarrow \int \mathrm{d} k.$

We define a operator ${a(k) = \frac{a_\lambda}{\sqrt{\Delta k}}}$ with ${k \approx k_\lambda}$ . The eletic field operator becomes

$\displaystyle \begin{aligned} E = & \sum_\lambda \sqrt{\Delta k} \sqrt{\frac{\hbar \omega_\lambda} {2\epsilon_0 LA}} \big( a(k_\lambda) e^{ik_\lambda z} + a^\dagger (k_\lambda) e^{-ik_\lambda z} \big) \\ =& \sum_\lambda \Delta k \sqrt{\frac{\hbar \omega_\lambda}{4\pi \epsilon_0 A }} \big( a(k_\lambda) e^{ik_\lambda z} + a^\dagger (k_\lambda) e^{-i k_\lambda z} \big) \\ \rightarrow & \int \mathrm{d} k \sqrt{\frac{\hbar \omega_\lambda}{ 4\pi \epsilon_0 A}} \big( a(k) e^{ikz} + a^\dagger (k) e^{-ikz} \big). \end{aligned} \ \ \ \ \ (2)$

${a(k), a^\dagger(k)}$ are continuous, have units of square root meters. The commutation is

$\displaystyle [ a(k), a^\dagger (k')] = \delta (k-k').$

The field Hamiltonian is

$\displaystyle H = \int \mathrm{d} k \hbar ck (a^\dagger (k) a(k) + \frac{1}{2} ).$

Now let’s transform the electric field into the spatial description.

$\displaystyle \left\{ \begin{array}{ll} a(z,t) =& \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} a(k,t) e^{ikz} \mathrm{d} k \\ a^\dagger (z,t) =& \frac{1}{\sqrt{2\pi}} \int_{-\infty}^{\infty} a^\dagger (k,t) e^{-ikz} \mathrm{d} k \end{array} \right.$

where the commutation is ${[a(z,t), a^\dagger(z',t)] = \delta (z-z')}$ . The dimension of ${a(z,t)}$ and ${a^\dagger(z,t)}$ is ${\frac{1}{\text{Length}}}$ . ${a^\dagger(z,t) a(z,t) \mathrm{d}z}$ means the number of photons in the space between ${z}$ and ${z+\mathrm{d}z}$ .

1.2. Maxwell-Bloch equation

The operator ${a(z,t)}$ ‘s motion follows Heisenberg equation. The Hamiltonian is ${H= H_{\mathrm{R}} + H_{\mathrm{atom}} + H_{\mathrm{int}}}$ . ${H_{\mathrm{R}}}$ represents the Hamiltonian of pure radiation field. ${H_{\mathrm{atom}}}$ is the Hamiltonian of matter. ${H_{\mathrm{int}}}$ is the interaction Hamiltonian between light and matter. Therefore we get

$\displaystyle \begin{aligned} \frac{\partial }{\partial {t}} a(z,t) =& \frac{1}{\sqrt{2\pi}} \int \frac{\partial}{\partial t} a(k,t) e^{ikz} \mathrm{d}k \\ =& \frac{1}{\sqrt{2\pi}} \int \frac{1}{i\hbar} \big[ a(k,t), H_{\mathrm{R}} + H_{\mathrm{int}} \big] e^{ikz} \mathrm{d}k. \\ \big[ a(k,t), H_{\mathrm{R}}] = &\int \mathrm{d}k' \hbar c k' \big[ a(k,t), a^\dagger (k',t) \big] a(k',t) \\ =& \hbar ck a(k,t). \end{aligned} \ \ \ \ \ (3)$

As ${\frac{\partial}{\partial z} a(z,t) = \frac{1}{ \sqrt{2\pi}} \int ik a(k,t) e^{ikz} \mathrm{d}k}$ , we get

$\displaystyle \frac{1}{\sqrt{2\pi}} \int \big[ a(k,t), H_{\mathrm{R}} \big] e^{ikz} \mathrm{d}k = -c \frac{\partial}{\partial z} a(z,t). \ \ \ \ \ (4)$

Combining Eq. (3) and (4) we get

$\displaystyle \big( \frac{\partial}{\partial t} + c \frac{\partial}{ \partial z} \big) a (z,t) = \frac{1}{i\hbar} \big[ a(z,t),H_{\mathrm{int}} \big]. \ \ \ \ \ (5)$

This is what to say Maxwell-Bloch equation, which describes the light field affected by atoms through the interaction ${H_{\mathrm{int}}}$ .

If the radiation field is restricted to the narrow band, that is ${\omega \approx \omega_0}$ , we get

$\displaystyle \begin{aligned} E=& \int \mathrm{d}k \sqrt{\frac{\hbar \omega_0}{4\pi \epsilon_0 A}} \big( a(k,t) e^{ikz} + a^\dagger (k,t) e^{-ikz} \big) \\ =& \sqrt{\frac{\hbar \omega_0}{2\epsilon_0 A}} \big( a(z,t) + a^\dagger (z,t) \big). \end{aligned} \ \ \ \ \ (6)$

In free space, we get time and space dependence ${a(z,t)= a(0,t-\frac{z}{c})}$ . Let’s define ${a(t)= \sqrt{c} a(z,t)}$ , ${a^\dagger(t) a(t)}$ is the flux of photons at time ${t}$ . For the light in the vacuum state we get

$\displaystyle \big[ a(t), a(t') \big] = \delta (t-t').$

1.3. Continuous Matter Operators

For simplicity, we only consider the collection of two-level atoms, which are coupled to the light by dipole transition. ${|g\rangle}$ and ${|e \rangle}$ are ground and excited states. The interaction Hamiltonian is

$\displaystyle H_{\mathrm{int}} = - \sum_j \vec{d}_j \cdot \vec{E} (\vec{R}_j), \ \ \ \ \ (7)$

where ${\vec{d}_j = -e \vec{r}_j}$ is the dipole operator for the ${j}$ ‘th atom, ${\vec{R}_j}$ is the location of the ${j}$ ‘th atom.

If we consider the light is linearly polarized along ${x}$ , we get Hamiltonian

$\displaystyle H_{\mathrm{int}} = \sqrt{\frac{\hbar \omega_0 }{2\epsilon_0 A}} \sum_j \big( d^* \sigma_{eg}^{(j)} a(z_j,t) + d \sigma_{ge}^{(j)} a^\dagger (z_j,t) \big). \ \ \ \ \ (8)$

It is convienient to define continuous spin operator

$\displaystyle \sigma_{ge} (z,t)= \frac{1}{\rho A \mathrm{d}z} \sum_{z \in [z,z+\mathrm{d}z ]} \sigma_{ge}^{(j)} (t), \ \ \ \ \ (9)$

where ${\rho A \mathrm{d}z}$ is the number of atoms in the slice ${[z,z+\mathrm{d} z]}$ , ${\sigma_{ge}(z,t)}$ is dimensionless. Rewriting Hamiltonian via continuous spin operator, we get

$\displaystyle H_{\mathrm{int}} = \sqrt{\frac{\hbar \omega_0}{2\epsilon_0 A}} \int_0^L \big( d^* \sigma_{eg} (z,t) a(z,t) + d \sigma_{ge} (z,t) a^\dagger (z,t) \rho A \mathrm{d}z. \ \ \ \ \ (10)$

Combining Maxwell-Bloch equation, we get

$\displaystyle \big( \frac{\partial}{\partial t} + c \frac{\partial}{\partial z} \big) a(z,t) = -i g \rho A \sigma_{ge} (z,t), \ \ \ \ \ (11)$

where ${g= \sqrt{\omega_0/2\omega_0\hbar A} d}$ . Together with the Heisenberg equation of motion for ${\sigma_{ge}(z,t)}$ , we have the coupled quantum Maxwell-Bloch equations describing light/matter interaction.

We have following commutation

$\displaystyle [\sigma_{\mu\nu}(z,t), \sigma_{\mu'\nu'}(z',t)] = \frac{1}{\rho A} ( \sigma_{\mu\nu'} \delta_{\nu\mu'} - \sigma_{\mu'\nu} \delta_{\nu'\mu} ) \delta(z-z').$

The dimension of ${\delta(z-z')}$ is ${\frac{1}{\text{Length}}}$ .

2. QND measurement induced entanglement

In this section we will discuss the scheme to realize QND measurement on matter and using the measurement to generate entanglement between atomic ensembles. The results we get in the last section will be widely used in this section.

The level struction of the atoms and schematic setup is shown in Lu-Ming Duan et al. PRL 85, 5643(2000).

The atomic ensemble is of a pencil shape with Fresnel number ${F= A/\lambda_0 L \approx 1}$ . The input laser pulse is linearly polarized with to orthogonal mode.

$\displaystyle E^+ = (z,t)= \sqrt{\frac{\hbar \omega}{4\pi \epsilon_0 A}} \sum_{i=1,2} a_i (z,t) e^{i(k_0 z - \omega_0 t)}, \ \ \ \ \ (12)$

where ${\omega_0 = \frac{2\pi c}{\lambda_0}}$ is the frequency of the laser. ${a^\dagger_i (z,t) a_i (z,t) \mathrm{d}z}$ is the photon number of ${i}$ th mode in the area ${[z, z+\mathrm{d}z]}$ . The dimension of ${a(z,t)}$ is ${\frac{1}{\sqrt{Length}}}$ . Suppose that ${\langle a_i(0,t) \rangle = \alpha_t}$ . The number of photons is ${2N_p = 2c \int_0^T | \alpha_t |^2 \mathrm{d}t \gg 1}$ .

Define Stokes operators

$\displaystyle \begin{aligned} S_x^p=& \frac{c}{2} \int_0^L (a_1^\dagger a_2 - a_2^\dagger a_1) d \tau, \\ S_y^p = & \frac{c}{2i} \int_0^L ( a_1^\dagger a_2 - a_2^\dagger a_1) d\tau, \\ S_z^p = & \frac{c}{2} \int_0^L (a_1^\dagger a_1- a_2^\dagger a_2) d \tau. \end{aligned} \ \ \ \ \ (13)$

The commutation of them is

$\displaystyle \begin{aligned} {S_x^P,S_y^p} =& \frac{ic^2}{2} \int_0^L \int_0^L \big[ a_1^\dagger a_1 \delta(\tau-\tau') - a_2^\dagger a_2 \delta(\tau-\tau') \big] d\tau d\tau' \\ =& \frac{ic^2}{2} \int_0^T (a_1^\dagger a_1- a_2^\dagger a_2) d\tau \\ =& i S_z^p. \end{aligned} \ \ \ \ \ (14)$

${[S^p_y, S^p_z] = iS_x^p}$ ; ${[S_z^p, S_x^p]= iS_y^p}$ .

For ouer coherent input, ${\langle S_x^p \rangle = N_p}$ , ${\langle S_y^p \rangle = \langle S_z^p \rangle = 0}$ . Because ${N_p \gg 1}$ , we can treat ${S_x^p}$ classically. We define ${X^p = \frac{S_y^p}{\sqrt{\langle S_x^p \rangle}}}$ , ${P^p = \frac{S_z^p}{\sqrt{\langle S_x^p \rangle}}}$ . We get ${[X^p, Y^p]=i}$ .

We use continuous atomic operators,

$\displaystyle \sigma_{\nu\mu} (z,t) = \lim_{\delta z \rightarrow \infty} \frac{1}{\rho A \delta z} \sum_{z \in (z,z+\delta z)} \sigma^{(j)}_{\nu\mu} (t), \ \ \ \ \ (15)$

where ${\nu,\mu = 1,2,3,4}$ . The Hamiltonian is

$\displaystyle H_I = \hbar \sum_{i=1,2} \int_0^L \big[ g \sigma_{i,i+2} a_i (z,t) + \mathrm{H.C.} \big] \rho A \mathrm{d} z, \ \ \ \ \ (16)$

where ${g= \sqrt{\frac{\omega_0}{2\omega_0 A \hbar}} d}$ . In the large detuing limit, we can adiabatically eliminated the upper state population from the equation of motion.

$\displaystyle \begin{aligned} H_{eff} =& \sum_{i=1,2} \hbar \int_0^L \big[ \frac{g^2}{\Delta} \sigma_{ii} a^\dagger_i (z,t) a_i (z,t) \rho A \mathrm{d}z \big]\\ =& \int_0^L \frac{g^2}{\Delta} \rho A \big[ (\sigma_{11} - \sigma_{22}) (a_1^\dagger a_1 - a_2^\dagger a_2) + \sigma_{11} a_2^\dagger a_2 + \sigma_{22} a_1^\dagger a_1 \big] \mathrm{d} z \\ =& \frac{g^2 \rho A}{\Delta} \int_0^L \big[ \frac{1}{2} (\sigma_{11} - \sigma_{22} )( a^\dagger_1 a_1 - a_2^\dagger a_2) + \frac{1}{2} N_a N_p \big] \mathrm{d} z, \end{aligned} \ \ \ \ \ (17)$

where ${N_a}$ is the number of atoms. Neglect the constant energy part, which only induces the global phase, The effective Hamiltonian is

$\displaystyle H_{eff} = \hbar \frac{2g^2}{\Delta L} \sqrt{N_aN_p} P^a P^p. \ \ \ \ \ (18)$

The Heisenberg equation is

$\displaystyle \begin{aligned} \dot{X}^p = &\frac{1}{i\hbar} \big[X^p, H_{eff}] \\ =& \frac{2g^2}{\Delta L} \sqrt{N_a N_p} P^p. \end{aligned} \ \ \ \ \ (19)$

We can eayily get the input-output relation as

$\displaystyle \begin{aligned} {X^P}' =& X^P + \frac{2g^2}{\Delta L} \sqrt{N_aN_p} P^a \frac{L}{c} \\ =& X^P + \frac{2g^2}{\Delta c} \sqrt{N_aN_p} P^a, \\ {X^a}' =& X^a + \frac{2g^2}{\Delta c} \sqrt{N_aN_p} P^p, \\ {P^\beta}' =& P^\beta, \quad \beta=a,p. \end{aligned} \ \ \ \ \ (20)$

Here we neglect the loss terms, which is small.

Using input-output relation (20), we get

$\displaystyle \begin{aligned} {X^P_2}' =& X_1^P + \kappa \big( P^a_1 + P^a_2 \big), \\ {X_1^a}' + {X_2^a}' = & X_1^a + X_2^a + 2\kappa P^p. \end{aligned} \ \ \ \ \ (21)$

where ${\kappa = \frac{2g^2}{\Delta c} \sqrt{N_aN_p}}$ . Is is easy to get the viarance of the operators

$\displaystyle \begin{aligned} \delta({X_2^P}') =& \frac{1}{2} + \kappa^2 \\ \delta({X_1^a}' + {X_2^a}') = &\delta(X_1^a) + \delta(X_2^a) + \delta(P^p) (2\kappa)^2 \\ =& \frac{1}{2} + \frac{1}{2} + 4\kappa^2 \frac{1}{2} \\ =& 1+ 2\kappa^2. \end{aligned} \ \ \ \ \ (22)$

Here we use the fluctuation ${\delta(X_1^a) = \delta(X_2^a) = \delta(P^p) =\frac{1}{2}}$ .

Suppose that ${\langle S_x^a \rangle_1 = \langle S_x^a \rangle_2 = N_a}$ , we get

$\displaystyle \big[ X_1^a + X_2^a, P_1^a + P_2^a ] = 2i. \ \ \ \ \ (23)$

$\displaystyle \Rightarrow \delta(X_1^a + X_2^a) \delta(P_1^a + P_2^a) =1.$

After measurement on ${{X^p}'}$ , the pile up noise in the variable ${X_1^a + X_2^a}$ is ${(1+2\kappa^2)}$ . The uncertainty relation must be hold on. Therefore, ${\delta(P_1^a + P_2^a) = \frac{1}{1+2\kappa^2}}$ . We rotate the collective atomic spins, ${X_1^p \rightarrow -P_1^a, P_1^a \rightarrow X_1^a , X_2^a \rightarrow P_2^a, P_2^a \rightarrow -X_2^a}$ . Then we measure again. The QND measurement will be on ${X_1^a - X_2^a}$ and makes it squeezed, ${\delta(X_1^a - X_2^a) = \frac{1}{1 + 2\kappa^2}}$ . We have prepared a two-mode squeezed state. The squeezing parameter ${r}$ is given by

$\displaystyle r = \frac{1}{2} \ln (1+ 2\kappa^2)$

Bose-Hubbard model in CQED systems

In this note, I will show how to realize effective Bose-Hubbard model in CQED systems. This note is mainly based on the paper Nature Physics 2, 849 – 855 (2006).

1. Introduction

Perturbative method can be well described the weak interation systems. But it is not useful in strong correlated systems. Besides, some key phenomenas in strong correlated systems are hardly ovserved because the time scale associated with them is very short. This is also the obstacles for realizing quantum processes.

Several important processes have been archieved in cold atomic systems, such as creating effective Bose-Hubbard model and realizing quantum phase transition in optical lattice trapped cold atoms. In this systems, a precise control of collective parameters is feasible. But it is still hard to indiviually address and control on single sites. Here a different approach to creat effective Bose-Hubbard model is proposed.

We consider an array of cavities. The dynamics of polaritons, combined atom photon excitations, in this setup is studied. Since the distance between adjucent cavities is considerably larger than the optial wavelength of the resonant mode., individual cavities can be addressed.

Model is shown bellow (refer to Fig 1 of quant-ph/0606097)

Photon hopping occurs between neighboring cavities. The repulsive force betweeen two polaritons occupying the same site is generated by a large Kerr nonlinearity. Three species of polaritons appear in our system. The species which is most steable is gonverned by effective Bose-Hubbard Hamiltonian.

$\displaystyle H_{eff} = \kappa \sum_{\vec{R}} (P_{\vec{R}}^\dagger)^2 (P_{\vec{R}})^2 + J \sum_{<\vec{R},\vec{R'}>} (P^\dagger_{\vec{R}} P_{\vec{R}'} + \mathrm{H.C.}), \ \ \ \ \ (1)$

where ${P^\dagger_{\vec{R}}}$ creates a polariton in the cavity of site ${\vec{R}}$ , and the parameters ${\kappa}$ and ${J}$ describe on site repulsion and inter cavity hopping respectively.

2. Quantised field in a periodic array of cavities

In periodic array of cavities, we have periodic dielectric constant

$\displaystyle \epsilon(\vec{r}) = \epsilon(\vec{r}+ \vec{R}), \vec{R} = \vec{n} R, \ \ \ \ \ (2)$

where ${\epsilon}$ is real. Asorption is neglected. The field is represented by vector potential ${\vec{A}}$ and a scalar potential ${\Phi}$ . They obey gauge condition: ${\Phi=0}$ and ${\Delta\cdot (\epsilon (\vec{r}) \vec{A}) = 0}$ . Maxwell equation is

$\displaystyle \begin{aligned} \vec{B}= \nabla \times \vec{A}; \quad \vec{E} = - \nabla \Phi - \frac{1}{c} \frac{\partial \vec{A}}{\partial t};\\ \vec{B}=\vec{H}; \quad \frac{1}{c} \frac{\partial \vec{D}}{\partial t} = \nabla \times \vec{H}. \end{aligned} \ \ \ \ \ (3)$

Let’s set ${\vec{D}= \epsilon(\vec{r}) \vec{E}}$ . Because of gauge choice, we get ${\vec{E}= -\frac{1}{c} \frac{\partial \vec{A}}{\partial t}}$ . Then we get

$\displaystyle \frac{\epsilon(\vec{r})}{c^2} \frac{\partial^2 \vec{A}}{\partial t^2} + \nabla \times (\nabla \times \vec{A}). \ \ \ \ \ (4)$

Now we discuss under which condition the tight-binding limit can be resonable. The design of the cavities must fulfilled the following three conditions. {}{}

(1) The resonant frequency of the caivities must falls within the forbidden gap of the surrounding struction. The coupling between cavities is due to the evanescent Bloch waves.

(2) Cavities are weakly coupled as the cavities are sufficently separated.

(3) The eigenmode of the field in a cavity of the system will remain nearly the same as the mode in a single cavity. Therefore, only nearest neighbor coupling between cavity modes must be taken into account.

Using tight-binding approximation, ${\vec{A}}$ can be expanded in Wannier function ${\omega_\Omega (\vec{r},t)}$ , each located at one single cavity at location ${\vec{R}}$ . Besides, it can also be expanded in Bloch wave ${\vec{A}_{\vec{k}} (\vec{r},t)}$ . The relation between Wannier function and Bloch function is

$\displaystyle \vec{A}_{\vec{k}} (\vec{r},t) = A_0 e^{i\omega_{\vec{k}} t} \sum_n e^{-ink\vec{R}} \vec{\omega}_{\Omega} (\vec{r} - nR\vec{e}_z). \ \ \ \ \ (5)$

From Eq. (3), we get

$\displaystyle \begin{aligned} \nabla \times (\nabla \times \vec{A}_{\vec{k}} (\vec{r},t) ) &= \epsilon( \vec{r}) \frac{\omega_{\vec{k} }^2}{c^2} \vec{A}_{\vec{k}} (\vec{r},t), \\ \nabla \times (\nabla \times \vec{\omega}_\Omega (\vec{r},t) )&= \epsilon_0 (\vec{r},t) \frac{\omega_c^2}{c^2} \vec{\omega}_\Omega (\vec{r},t). \end{aligned} \ \ \ \ \ (6)$

${\omega_{\Omega} (\vec{r},t)}$ is real and normalized, ${ \int \mathrm{d}^3 \vec{r} \epsilon_0 (\vec{r}) \vec{\omega}_{\Omega} (\vec{r}) \cdot \vec{\omega}_{\Omega} (\vec{r}) = 1}$ . Multiplying both sides of Eq. (6) from left-hand side by Eq. (5) and spatailly integrating, we get the left hand-side of Eq. (6) becomes

$\displaystyle \begin{aligned} &\int \mathrm{d}^3 \vec{r} \epsilon(\vec{r}) \frac{ \omega_{ \vec{k}}^2}{c^2} A_0 e^{i\omega_{\vec{k}}t} \sum_n e^{in K \vec{R}} \vec{\omega}_\Omega (\vec{r}) \cdot \vec{\omega}_\Omega (\vec{r} - nR\vec{e}_z) \\ =&\big[ 1 + \Delta \alpha \sum_{n\neq 0} e^{-inkR} \alpha_n \big] \frac{\omega_{k}^2 }{c^2} A_0 e^{i\omega_k t}, \end{aligned} \ \ \ \ \ (7)$

where ${\alpha_n = \int \mathrm{d}^3 \vec{r} \epsilon(\vec{r}) \vec{\omega}_{\Omega} (\vec{r}) \cdot \vec{\omega}_\Omega (\vec{r} - n R\vec{e}_z)}$ and ${\Delta \alpha = \int \mathrm{d}^3 \vec{r} [\epsilon(\vec{r}) - \epsilon_0 (\vec{r}) ] \vec{\omega}_\Omega (\vec{r}) \cdot \vec{\omega}_\Omega (\vec{r})}$ . The right hand-side of Eq. (6) reads

$\displaystyle \begin{aligned} &\int \mathrm{d}^3 \vec{r} \vec{\omega}_{\Omega} (\vec{r}) \cdot \nabla \times (\nabla \times \vec{A}_{\vec{k}} )\\ =&\big[ 1+ \sum_{n\neq 0} e^{-in kR} \beta_n \big] A_0 e^{i\omega_{\vec{k}} t} \frac{\omega_c^2}{c^2} , \end{aligned} \ \ \ \ \ (8)$

where ${\beta_n = \int \mathrm{d}^3 \vec{r} \epsilon_0 (\vec{r}- nR \vec{e}_z )\vec{\omega}_\Omega (\vec{r}) \cdot \vec{\omega}_\Omega (\vec{r}- nR \vec{e}_z)}$ . Combining Eq. (7) and (8), we get dispersion relation

$\displaystyle \omega_{\vec{k}} \big[ 1+ \Delta \alpha + \sum_{n\neq 0} e^{-in kR} \alpha_n \big] = \omega_c^2 \big[ 1+ \sum_{n\neq 0} e^{-inkR} \beta_n \big]. \ \ \ \ \ (9)$

If the coupling between the resonantor is sufficiently weak, we can count only the nearest neighbor coupling, ie. ${\alpha_n = \beta_n = 0}$ , if ${n\neq \pm 1}$ . For symmetry considerations, we also require that ${\alpha_1 = \alpha_{-1}}$ , ${ \beta_1 = \beta_{-1}}$ . Finally we assume ${\alpha_1}$ , ${\beta_1}$ and ${\Delta \alpha}$ to be small, we can simplify the dispersion relation to

$\displaystyle \omega_{k} = \omega_c \big[ 1 - \frac{\Delta \alpha}{2} + \kappa_1 \cos (kR) \big], \ \ \ \ \ (10)$

where we define the coupling factor ${\kappa_1}$ as ${\kappa_1 = \beta_1 - \alpha_1}$ . Suppose that ${\Delta \alpha}$ is neligiable compared with ${\kappa_1}$ . The dispersion relation can be simplified as

$\displaystyle \omega_k = \omega_c [ 1+ \kappa_1 \cos (kR) ]. \ \ \ \ \ (11)$

Now following the usual way to quantize the system from the eigenmodes ${\omega_k}$ , we can get Hamiltonian of the system

$\displaystyle H = \omega_c \sum_k b^\dagger_k b_k [ 1+ \kappa_1 \cos (kR)]. \ \ \ \ \ (12)$

From now on, we use natural unit that ${\hbar = 1}$ . ${b^\dagger_k}$ and ${b_k}$ is the creation and annihilation operators of a photon occupying the extended eigenmode. It would be conviend to transform the operators in Wannier representation, where the operators describe the localized eigenmodes (cavity modes). The transformation connections is

$\displaystyle b_k = \frac{1}{\sqrt{N}} \sum_n e^{-inkR} a_n, \ \ \ \ \ (13)$

where ${a_n= a(nR)}$ . In the localized modes, the Hamiltonian (12) reads

$\displaystyle \begin{aligned} H= &\frac{\omega_c}{N} \sum_k \sum_n \sum_{n'} e^{ikR(n-n')} a^\dagger_n a_{n'} \big[ 1 + \kappa_1 \cos(kR)\big]\\ =& \omega_c \sum_{n=1}^N a_n^\dagger a_n + \frac{\kappa_1 \omega_c }{2N} \sum_{n,n'} \sum_{k} e^{ikR(n-n')} a^\dagger_n a_n' (e^{ikR}+ e^{-ikR}) \\ =& \omega_c \sum_{n=1}^N a_n^\dagger a_n + \frac{\kappa_1 \omega_c}{2N} \sum_{n,n'} \sum_k \Big[ e^{ikR(n-n'+1)} a^\dagger_n a_n' + e^{ikR(n-n'-1)} a^\dagger_n a_n' \Big] \\ =& \omega_c \sum_{n=1}^N a_n^\dagger a_n + \frac{\kappa_1 \omega_c}{2} \sum_{n=1}^N [a^\dagger_n a_{n+1} + \mathrm{H.C.} ] + \frac{N\kappa_1 \omega_c}{2} \end{aligned} \ \ \ \ \ (14)$

Neglect the total energy term ${\frac{N\kappa_1 \omega_c}{2}}$ , we get the effective Hamiltonian is

$\displaystyle H_{eff} = \omega_c \sum_{n=1}^N a_n^\dagger a_n + \frac{\kappa_1 \omega_c}{2} \sum_{n=1}^N [a^\dagger_n a_{n+1} + \mathrm{H.C.} ], \ \ \ \ \ (15)$

where we use the donate ${a_{N+1} = a_N}$ .

3. Site repulsion

To generate a repulsion between polaritons that are located in the same cavity, we fill the cavity with 4 level atoms.

In a rotating frame with ${H_0 = \omega_c (a^\dagger a + 1/2) + \sum_{j=1}^N (\omega_c \sigma_{22}^j + \omega_{c} \sigma_{33}^j + 2 \omega_c \sigma_{44}^j)}$ , the Hamiltonian of the atoms in the cavity reads

$\displaystyle H_I = \sum_{j=1}^N \big[\varepsilon \sigma_{22}^j + \delta \sigma_{33}^j + (\Delta + \varepsilon) \sigma_{44}^j \big] + \sum_{j=1}^N ( \Omega_L \omega_{23}^j + g_{13} \sigma_{13}^j a^\dagger + g_{24} \sigma_{24}^j a^\dagger + \mathrm{H.C.} ), \ \ \ \ \ (16)$

where ${\sigma_{kl}^j = |k_j\rangle \langle l_j|}$ . ${\omega_c}$ is the frequency of the cavity mode. ${\Omega_L}$ is the Rabi frequency of the driving field. ${g_{13},g_{24}}$ are the parameters of the dipole coupling of the cavity mode.

The Hamiltonian truncated to the subspace at most two excitations can be diagonalized. Let’s define the following operators ${ S^\dagger_{12} = \frac{1}{\sqrt{N}} \sum_{j=1}^N \sigma_{21}^j}$ and ${ S^\dagger_{13} = \frac{1}{\sqrt{N}} \sum_{j=1}^N \sigma_{31}^j}$ . In the limit ${N\gg 1}$ and excitation is restricted to ${2}$ , we get $latex \displaystyle \begin{aligned} \big[ S_{12}, S^\dagger_{12} \big] =& \frac{1}{N} \sum_{j=1}^N |1\rangle \langle 1|_j \approx 1 ,\\
\big[ S_{13}, S^\dagger_{13}\big] = & \frac{1}{N} \sum_{j=1}^N |1\rangle \langle 1|_j \approx 1 ,\\
\big[ S_{12}, S^\dagger_{12} \big] = & \big[ S_{13}, S^\dagger_{13} \big] =0 \\
\big[S_{12}^\dagger, S_{13} \big] = &\big[S_{13}^\dagger, S_{13} \big] = 0. \end{aligned}$ Therefore, ${S_{12}}$ and ${S_{13}}$ fulfill boson commutation and represent two different bosons.

Neglecting two photon detuning and coupling to level ${4}$ . Setting ${g_{24}=0}$ and ${\varepsilon =0}$ , level ${4}$ decouples from the rest of the two exciation manifold. We assume the number of atoms is large, ${N \gg 1}$ . The Hamiltonian (16) in terms of ${S_{12}}$ and ${S_{13}}$ reads

$\displaystyle H_I = \delta S^\dagger_{13} S_{13} + \Omega_L S^\dagger_{12} S_{13} + g S_{13} a^{\dagger} + \Omega_L S_{13}^\dagger S_{12} + g S_{13}^\dagger a, \ \ \ \ \ (17)$

where ${g= \sqrt{N} g_{13}}$ . Let’s define operator ${Q_0^\dagger = \frac{1}{B} ( \Omega_L S_{12}^\dagger + ga^\dagger )}$ , where ${B= \sqrt{\Omega_L + g^2}}$ . ${Q_0}$ is boson which fulfills ${[Q_0, Q^\dagger_0] = 1}$ . We can define another boson operator ${P_0^\dagger= \frac{1}{B} ( g S_{12}^\dagger - \Omega_L a)}$ which fulfils ${[P_0,P^\dagger_0]=1}$ , ${[P_0, Q_0] = [P_0, Q_0]=0}$ . So ${P_0}$ and ${Q_0}$ are two different bosons.

In terms of ${Q_0}$ and ${S_13}$ , Hamiltonian (16) reads

$\displaystyle H_I = \delta S_{13}^\dagger S_{13} + B (Q^\dagger_0 S_{13}+ S^\dagger_{13} Q_0). \ \ \ \ \ (18)$

We want use normal transformation to diagonlize ${H_I}$ . Let’s define

$\displaystyle \begin{aligned} P^\dagger_+ = &\mu S^\dagger_{13} + \nu Q^\dagger_0 \\ P^\dagger_+ =& -\nu S^\dagger_{13} + \mu Q^\dagger_0, \end{aligned} \ \ \ \ \ (19)$

where ${\nu}$ and ${\mu}$ are parameters which fulfill ${|\nu|^2 + |\mu|^2 = 1}$ . Is is easy to prove that ${[P_+ ,P^\dagger_+] = [P_-, P^\dagger_-] = 1}$ and ${[P_+, P_-^\dagger] = [P_-, p_+^\dagger] =0}$ . So ${P_-}$ and ${P_+}$ are two different bosons. ${H_I}$ in terms of ${P^\dagger_+}$ and ${P^\dagger_-}$ reads

$\displaystyle \begin{aligned} H_I = &\delta (\mu P^\dagger_+ - \nu P^\dagger_-) (\mu P_+ - \nu P_-) + B\big[ (\nu P^\dagger_+ + \mu P^\dagger_-) ( \mu P_+ - \nu P_-) \\&+ (\nu P^\dagger_+ - \nu P^\dagger_-) ( \nu P_+ + \mu P_-) \big] \end{aligned} \ \ \ \ \ (20)$

To diagonlize ${H_I}$ , ${\nu}$ and ${\mu}$ must fulfill

$\displaystyle B(\nu^2 - \mu^2) + \delta \mu\nu = 0.$

Solve equation we get ${\frac{\nu}{\mu} = \frac{\delta \pm A}{2B}}$ , where ${A = \sqrt{\delta^2 + 4B^2}}$ . The diagonlized Hamiltonian is

$\displaystyle H = \mu_0 P^\dagger_0 P_0 + \mu_+ P^\dagger_+ P_+ + \mu_- P^\dagger_- P_-, \ \ \ \ \ (21)$

where ${\mu_0=0}$ , ${\mu_+ = (\delta+A)/2}$ and ${\mu_- = (\delta-A)/2}$ .

To write the full Hamiltonian ${H_I}$ in the polariton basis, we express operators ${\sum_{j=1}^N \sigma_{22}^j}$ and ${a^\dagger \sum_{j=1}^N \sigma_{24}^j}$ in terms ${P^\dagger_0}$ , ${P_+^\dagger}$ and

${P^\dagger_-}$ . Let’s suppose that ${|g_{24}|,|\varepsilon|, |\Delta| \ll |\mu_+ |, |\mu_-|}$ . Therefore rotating wave approximation is used here. The interaction between ${P_0}$ and ${P_+}$ , ${P_-}$ can be negelected. Therefore we have

$\displaystyle \begin{aligned} \varepsilon \sum_{j=1}^N \sigma_{22}^j = &\varepsilon S_{12}^\dagger S_{12} \\ \approx & \varepsilon \frac{g^2}{B^2} P^\dagger P_0 \\ g_{24} \sum_{j=1}^N (\sigma_{42} a + \mathrm{H.C.} =& g_{24} (S^\dagger_{14} S_{12} a + \mathrm{H.C.} \\ \approx & - g_{24} \frac{g \Omega_L}{B^2} (S^\dagger_{14} P_0^2 + \mathrm{H.C.}). \end{aligned} \ \ \ \ \ (22)$

At last, level ${4}$ can be adiabatically eliminated. The Hamiltonian we get at last is

$\displaystyle H = \kappa (P_0^\dagger)^2 P^2, \ \ \ \ \ (23)$

where

$\displaystyle \kappa = \frac{g_{24}^2}{\Delta} \frac{N g_{13} \Omega_L}{(N g_{13}^2 + \Omega_L^2 )^2}$

. For simplicity, we set ${\varepsilon =0}$ .

4. The total Hamiltonian

If we use polariton operators to write the hopping terms in Eq. (1),

$\displaystyle a^\dagger_n a_{n+1} \approx \frac{\Omega_L^2}{B^2} P^\dagger_n P_{n+1}.$

Contribution of different polaritons decouple due to the separation of their energy. There for we get the Hamiltonian ${P_n^\dagger}$ takes the on the form (1), with

$\displaystyle J= \frac{ \omega_C \Omega^2_L}{2(N g_{13}^2 + \Omega_L^2)} \kappa_1.$

一维薛定谔方程的数值解法

1. 物理想法

参照《计算物理学》(Steven E.Koonin著，秦克诚译)第三章的相关内容，我写了这篇笔记。

我们考虑的问题是求一个在一维位势~ ${V(x)}$ ~中运动的质量为~ ${m}$ ~的粒子的量子力学定态解。我们假定~ ${V(x)}$ ~的形状如图，是一个三角形势阱，

在~ ${x=x_{min}}$ ~与~ ${x=x_{max}}$ ~两点处位势变为无穷，这两点之间则有一个势阱。我们假定定态不含时的薛定谔方程和边界条件是

$\displaystyle \begin{aligned} \frac{\mathrm{d} \phi}{\mathrm{d} x^2} + k^2(x) \phi(x) = 0; \\ \psi(x_{min}) = \psi(x_{max})=0, \end{aligned} \ \ \ \ \ (1)$

其中

$\displaystyle k^2(x)= \frac{2m}{\hbar^2} [E - V(x)].$

我们要求是非零的能量本征值。对这个本征值~ ${E}$ ~来说，在经典力学容许的范围内( ${E<V(x)}$ )是振荡的，在经典力学禁戒的区域( ${E>V(x)}$ )呈现指数行为。整体来看，本征值小于零的解是束缚态的解，本征值大于零的解是连续解。我们这里感兴趣的是束缚态。

要算出这个方程束缚态的本征值，我们需要利用打靶法计算。由于积分区域既有振荡部分，也有指数部分，直接积分会引起数值上的不稳定。所以我们不能把波函数直接从区域的左边积到右边，再利用边界条件判断是否是本征值。可是考虑到波函数的连续性以及其导数的连续性，我们可以分别从左边积分到右边，再从右边积分到左边，在振荡区和指数区的交界处，我们来判断波函数及其导数是否连续。考虑到波函数归一化时我们总可以选择一种归一化系数，使得两种不同的积分方式得到的波函数大小在交界处的值相等，所以我们只需要保证导数连续即可。也就是下式要成立：

$\displaystyle \frac{\mathrm{d}\psi_<}{\mathrm{d}x}\Bigg\vert_{x_m}- \frac{\mathrm{d}\psi_>}{\mathrm{d}x}\Bigg\vert_{x_m} = 0. \ \ \ \ \ (2)$

这里的~ ${\psi_<}$ ~和~ ${\psi_>}$ ~分别是向右和向左积分得到的波函数。用最简单的有限差分格式代替这个导数式，我们得到下式：

$\displaystyle f \equiv \frac{1}{\psi} [\psi_<(x_m - h) - \psi_> (x_m - h)] = 0, \ \ \ \ \ (3)$

由于归一化保证了~ ${\psi_< (x_m) = \psi_> (x_m)}$ ，上式中的~ ${\psi}$ ~可以用~ ${x_m}$ ~点上的~ ${\psi_<}$ ~来代替。如果不存在交界区(如具有上图所示位势，而~ ${E>0}$ )，那么~ ${x_m}$ ~可以选择在任何区域，如果有多余两个转折点，那么就必须把每个区域中的解拼在一起。

之所以用这种拼凑的解法，是因为我们只能保证积分在每个小区域中是准确的，但是多个区域之间利用连续性条件我们可以凑成一个完整的波函数。

2. 算法分析

我们为了获得精确的积分结果，用到了~Numerov~算法。对薛定谔方程来说，我们利用~Numerov~算法得到递推关系如下：

$\displaystyle \Bigg(1+\frac{\hbar^2}{12} k^2_{n+1} \Bigg) \psi_{n+1} - 2\Bigg( 1 - \frac{5\hbar^2}{12} k^2_n \Bigg) \psi_n + \Bigg( 1 + \frac{\hbar^2}{12} \Bigg) \psi_{n-1} = O(h^6). \ \ \ \ \ (4)$

为了计算积分，我们必须利用泰勒级数算出头两项的数值出来。由于波函数的归一是任意的，所以波函数在边界处的一阶导数可以随我们任意取，不影响最后得到的本征值和波函数结果。

为了让计算机更好处理，我们把薛定谔方程无量纲化，写成

$\displaystyle \Bigg[ - \frac{1}{\gamma^2} \frac{\mathrm{d}^2}{\mathrm{d} x^2} + v(x) - \epsilon \Bigg] \psi(x) = 0, \ \ \ \ \ (5)$

其中

$\displaystyle \gamma^2 = \frac{ 2m a^2 V_0}{\hbar^2}$

是系统经典力学本性的一个无量纲测度，~ ${\epsilon = E/V_0}$ ~是无量纲能量，坐标由长度~ ${a}$ ~标度。

我们用打靶法获得本征值。先从能量的最小可能取值开始往上搜索。如果发现~ ${f}$ ~变号，那么可以认为这里可能有一个本征态。于是粗略的获得本征值的范围，然后减小步长，获得更加精确的结果，直到~ ${f}$ ~的绝对值小于~ ${10^{-6}}$ ~为止。由于算法的问题，~ ${f}$ ~变号不能一定保证就有本征值在附近。为了排除这种情况，我们要联合考虑看是否随着搜索的进行~ ${f}$ ~的绝对值在不断减少。如果不是如此，那么我们就要去掉这个结果。在得到一个能量本征值之后，我们以这个能量本征值为新的起点，继续前面的过程，可以获得一系列的能量本征值~ ${E_i}$ ，直到能量本征值大于零为止。对每一个本征值对应的波函数，我们要分别归一化。归一化时，我们用~Simpson~法求积分，比梯形法具有更高的精确度。

我们考虑一个具体的问题：半导体硅中低能谷电子有效质量为~ ${m^*=0.916 m_e}$ ，势场宽为~ ${35 nm}$ ~,外电场从1伏特变到零，电子在外场中的势能从~ ${-1}$ ev~变到零。我们把积分区域分成500个格子，起步能量步长为~ ${1/500}$ ~电子伏特，能量尺度设为~0.1ev，长度标度为~ ${a=3}$ nm。最后算出能量本征值为~ -9.2426, -8.6758,-8.2117, -7.8016 , -7.4267, -7.0773 , -6.7477, -6.4340, -6.1336 ,-5.8444, -5.5649, -5.2940, -5.0307,-4.7741, -4.5238, -4.2790,-4.0393, -3.8044, -3.5738, -3.3473, -3.1246, -2.9055, -2.6897,-2.4770, -2.2673, -2.0604, -1.8562, -1.6544, -1.4551,-1.2581, -1.0633, -0.8705, -0.6794, -0.4889, -0.2965, -0.0984，这里能量的单位是~ ${0.1}$ ev。可以看出，随着能级的升高，能量之间的间隔越来越小。

为了确保算出的能量本征值是正确的，我们把两个不同能量本征值对应的波函数相乘再积分，发现二者是正交的，正好满足量子力学中不通能量本征值的波函数必须正交的性质。下面给出几个图，看看不同能量本征值下的波函数有什么特点，并与前面理论分析的结果相互比较，以确保我们的计算无误。

如图，左边是第三个能级的波函数，右边是第八个能级的波函数。可以看出能级数与振荡部分波峰波谷之和相等。在经典禁戒区，波函数是指数衰减的。与我们前面分析的结论是一致的。

下面附有源程序

clear; m0=9.108E-31; ml=0.916*m0; h1=1.055e-34; N=500; M=500; tol=1e-6; a=3e-9; V0=1.602e-20; l=35e-9; dx=l/N; L=l/a; Dx=dx/a; X=0:Dx:L; x=X*a; ga=sqrt(2*ml*V0*a^2/h1^2); Vm=1; V=-Vm:Vm/N:0; V=V*1.602e-19; v=V/V0; eold=v(1)+abs(v(1))/M; for k=1:N/2 U(k,2*k-1:2*k+1)=[1,4,1]; end for k=1:1000 e1=eold+abs(v(1))/M; err=2*tol; v1=-(v-eold); s=find(v1tol & abs(err)<10000 & e1<0 & M<10^10 v1=-(v-e1); s=find(v1<=0);

if length(s)==1 e=e(find(e~=0)); E=e*V0; P=Psi(find(e~=0),:); return end K=-ga^2*(v-e1); psi(N+1)=0; psi(N)=Dx*0.5; for n=N+1:-1:s(1)+1 psi(n-2)=(2*(1-K(n-1)*5*Dx^2/12)*psi(n-1)-(1+K(n)*Dx^2/12)… *psi(n))/(1+K(n-2)*Dx^2/12); end psi1(1)=0; psi1(2)=Dx*0.5; for n=1:s(1)-1 psi1(n+2)=(2*(1-K(n+1)*5*Dx^2/12)*psi1(n+1)-(1+K(n)*Dx^2/12)… *psi1(n))/(1+K(n+2)*Dx^2/12); end psi=psi/psi(s(1)); psi1=psi1/psi1(s(1)); err=psi1(s(1)+1)-psi(s(1)+1); eold=e1; e1=e1+abs(v(1))/M; if err*errold<0 eold=e1-3*abs(v(1))/M; M=M*10; e1=eold+abs(v(1))/M; end end if abs(err)<0.0001 e(k)=eold; end Psi(k,1:s(1))=psi1(1:s(1)); Psi(k,s(1)+1:N+1)=psi(s(1)+1:N+1); Psi(k,:)=Psi(k,:)/sqrt(Dx*sum(U*Psi(k,:).^2')/3); psi=0; psi1=0; M=500; eold=eold+2*abs(v(1))/M; end