我们生活在广漠宇宙的一隅，如果宇宙是片海洋，那么星系就是一个个岛屿。据估计，宇宙中有上千亿个星系（如果包含比较小的星系，最近的观测发现可能有上万亿个星系），而一个星系中又有上千亿颗恒星，由此看来，太阳只是我们银河系中普通的一颗恒星。当凝视一张星系的图片时，我们往往会被星系美丽结构所吸引。那么比星系更大的尺度上是否也存在结构呢？如果有，这些结构的形成又有怎样的历史呢？
一个世纪前，人类还不清楚自己的银河系大小，而现代的大型天文望远镜，可以把我们的目光带到宇宙的深处。图1有点像医学检查用的CT片，它来自于斯隆数字巡天（SDSS）项目，其中每一个点代表着一个星系。星系的分布之所以是两个扇面形状，是由于只观测到了这些区域，并不是扇形之外的区域没有星系存在。作为观测者，我们处于图片的中心，越远离中心的位置离我们也就更遥远。
在宇宙学中常用红移（redshift，一般记成z）来标记时间，红移越大，时间就越早，今天宇宙的红移为零。图1中的红移覆盖到0.15，这意味着本图看到了宇宙大约19亿年前发出的光。相比之下，银河系中心的光到达我们的眼睛大约需要3万年，而火星到地球的距离平均只有12个光分。
图中星系的分布明显是有结构的，总体来说像一张大网。
在有些地方聚集了大量的星系，有些地方则是星系的“空洞”，在星系的高度聚集区之间有很多丝状的结构。如何来理解宇宙在星系尺度之上的复杂结构呢？经典的理论分析能够在“线性”近似的层面上给出很多解释。所谓线性就是说等比例的改变初始条件，也会得到等比例改变的结果。然而，宇宙结构的形成是一个非线性的过程，而非线性的系统则不会遵守这样的“规则”，微小的初始的差别会造成后来巨大的不同。
因此，即使把理论计算改进得很复杂依然难以达到目前所需要的水平。由此，计算机数值模拟就成了一个帮助我们理解宇宙的不可或缺的工具。但是要开展数值模拟，我们必须首先提供宇宙的初始条件和动力学演化方程，下面我们逐一介绍。
那么宇宙结构的初始条件是什么样的呢？（注意我们讨论的是宇宙结构的开始，而不是宇宙本身的开始。）这必须要回到宇宙诞生的早期去。
根据暴胀理论，宇宙在经历过极早期的急速暴胀之后开始进入热大爆炸阶段，这时的宇宙已经被撑成一个极端平直的空间，同时尺度也远大于我们能够观测的宇宙视界的大小，接着物质开始在宇宙中形成。由于那时的宇宙能量极高，整个宇宙中的物质存于高温状态，物质之间频繁地相互作用。打个比方，宇宙就像是煮开的一锅粥，所以不难理解宇宙的结构在初始时是非常均匀的。
也就是说，那时的宇宙中很难有结构，即使产生了结构也会被“煮化了”。
关于这一点，宇宙微波背景辐射的观测提供了非常强的支持，宇宙微波背景辐射温度来自宇宙诞生之后37万年左右（红移大约在1100），它在各个方向上都极度均匀，从而得出：宇宙在大尺度上是均匀并且各向同性的，也就是说，宇宙中没有特殊的位置。当然自然界中不存在完美的平滑。宇宙的不均匀性在不同尺度上也是不一样的，尺度越小，不均匀的涨落程度越大。而对这个均匀背景的微小的偏离就能提炼出大量宇宙早期的信息。
由于宇宙在不断膨胀，物质能量密度也就越来越稀疏，最终炽热的原初火球逐渐冷却了下来，宇宙也进入了黑暗时期（Dark Ages）。伴随气体的凝聚，恒星和星系开始形成，宇宙再次被电离。没有了中性氢的阻挡，光子自由地在宇宙空间中穿行，又经过130亿年的时间，宇宙终于演化到了今天。
宇宙早期残留下的这些高高低低的密度涨落，就是未来宇宙大尺度物质结构的种子。也是宇宙学数值模拟的入口，通常我们就是从大爆炸之后的1~2千万年左右开始（红移z≈100）模拟的，那时宇宙的大小只有现在的百分之一左右。
明确了宇宙结构的种子，我们还需要了解控制结构演化行为的第二个因素——动力学。世界上有四种力，强相互作用和弱相互作用都是短程力，所以能在大尺度上起作用的只有电磁力和引力。通常电磁力要比引力强很多倍，但是宇宙不存在显著的净电荷或者净电流，所以电磁力在结构演化中也不重要，担当宇宙结构演化的主角的必然是引力。
按照广义相对论所言，引力在本质上是时空的弯曲，对于宇宙来说时空的弯曲可以理解成整体和局部两个部分。整体部分就是宇宙空间整体上在膨胀，这是有宇宙的全部物质能量共同决定的。具体来说就是，宇宙的物质能量包含暗能量、暗物质、重子物质和辐射四类组分。暗能量使得宇宙加速膨胀，我们虽然还不知道暗能量的本质，但是对于通常的模型在合理的假设下，它是无法形成明显的结构的。
辐射指的是类似光子或者中微子这样以接近光速运动的能量物质。它们速度太快，不会被束缚在局部引力场中，所以辐射也不会形成明显的结构。那么时空局部的弯曲实际上就决定于暗物质和重子的分布。原则上时空的弯曲需要依据广义相对论来计算，但是如此大规模的数值相对论的复杂度是现代计算机能力所无法做到承受的。
幸运的是，万有引力是广义相对论的弱场近似，而牛顿力学在太阳系内的应用已经常准确，因此，同样可以把它应用到比太阳系更加稀薄的宇宙空间（也就是更弱的引力场）。
宇宙密度的演化有两个互相竞争的过程：（1）宇宙膨胀带来的“减速”，想象一个不开推进器的飞船朝向目标前进，由于目标被宇宙膨胀拉得越来越远了，所以相对于这个参考点来说就是飞船的速度是“变慢了”；（2）物质间的引力，飞船会被物质密度高的地方（高于宇宙背景）吸引，反过来说就是被物质密度低于宇宙背景的地方“排斥”。
所以宇宙中的结构会在初始结构种子的基础上，在引力的作用之下，高密度的地方变得更加密集，低密度的地方变得更为稀疏（非常类似于人口流动：“有吸引力的地方变得更密集，而荒芜的地方变得更荒芜”），所以宇宙结构的演化是一个不断“增长”的过程。
根据这些分析，我们就可以开始模拟一个宇宙了。（1）由于光速的限制，我们所观测到的永远是这个宇宙（即便你相信多重宇宙存在）中的一部分，而且很可能是其中很小的部分。
所以我们可以取一个方盒子，其中盛入一部分宇宙。这个盒子必须足够大（比如几百Mpc以上，1 Mpc≈326万光年），这样才能代表宇宙的整体。如果盒子中的物质流出去怎么办？只要周期性地把流出的部分从盒子的另一边流回来，就能保证总的质量守恒。但是有一个问题，我们选取的这个盒子里的结构并不会碰巧与我们周围的宇宙一致（这个问题是内在的，原因在于初始条件信息的不完备）。
即便如此，这样的模拟依然是有意义的，因为宇宙学原理告诉我们，在不同的地方宇宙的统计性质是一样的。也就是说，我们虽然造不出一模一样的银河系和太阳，但是我们不难在盒子里找到一些相似的对应物。
（2）我们只要模拟物质就够了（因为暗能量和辐射都是不形成结构的），为了描述物质的分布，通常有两种方法，一种是把模拟的盒子分成网格，网格上计算能量、密度、速度场等等信息，演化遵守宇宙学框架下的流体力学；一种方法是把宇宙全部的质量分成一些粒子，每个粒子遵守膨胀空间的万有引力，也就是宇宙学N体模拟。后一种方法更为普遍一些。
这是因为我们关心的尺度通常从上千万光年到几千光年甚至更小，而计算的分辨率是按照最小的尺度来决定的（动态范围太大），而粒子位置的精度是很容易达到这个范围的（单精度的计算也能覆盖100万倍的动态范围），同时粒子数的守恒就天然保证总质量的守恒，这也是个很大的优势。
下面就来看一下数值模拟中结构形成的历史。
图4最上边是宇宙年龄14亿岁的时候（z=4.0，也就是今天宇宙尺度的1/5时），而最下边的是宇宙今天的样子，也就是年龄138亿左右（z=0）。从这六个时间切片可以看出，宇宙在初始密度稍微高的地方开始不断聚集，空旷的区域也更加明显。在密度最大的地方，物质形成了星系团这样巨大的结构（Mpc尺度），这个结构中物质的动能和引力势能平衡，密集程度不再明显。
但是它作为一个大势井，外部的物质还会不断地落入、合并。
接下来我们讨论一下丝状结构，这个结构在星系的观测（图1）和模拟中（图4）都是存在的。首先想象一个均匀对称的球体，在自引力的作用下发生塌缩，塌缩的速度在各个方向自然也是对称的，所以最终形成的结构大致也是球对称的。
现在考虑一个棒状系统的塌缩，因为加速度是按照距离平方减小，所以在棒的长轴方向的加速度大大弱于短轴，结果自然是短轴塌缩更快，而长轴方向塌缩较慢。宇宙结构的种子大致是一个高斯随机场，这意味着存在各种可能的结构，早先类似棒状的系统在今天看来就像是一条条的丝线。因此宇宙结构的图景就是像一张巨大的网络，在丝线连接的节点上是大的星系团。
N体模拟的原理非常简单，就是计算粒子对之间的引力，但是这在计算上的消耗是非常大的。最简单而又准确的方法就是累加所有粒子对之间的引力，然后根据这一时刻的力去更新粒子速度，进而用当前的速度去更新粒子的位置，在新的位置上重新计算新的力，如此循环往复。但使用这种直接法去计算所引力的计算时间复杂度是粒子数的平方O（N2），也就是计算100个粒子需要计算10000对，而计算1万个粒子则需要1亿对。
所以随着粒子数目增加，很快就会耗尽计算机的资源。
20世纪70年代皮布尔斯（Peebles）等就是用这样的方法进行了第一个宇宙学意义上的模拟，当时的粒子数只有1500个。随着计算机硬件的发展，cpu的运算能力越来越高，相同时间里计算的粒子数目相应地增加。即便如此，在80年代之前，使用直接计算引力对的方式的模拟所能使用总粒子数依然没有超过1万。真正让宇宙学模拟的粒子数提高所依靠的是算法的改进。
到了80年代，引入了粒子网格（Particle-Mesh）近似计算方法，粒子数目上升到万级。
但是这一方法对于粒子高度成团的时候是无效的，这时计算引力的树方法（Tree）和自适应网格加密（Adaptive Mesh Refinement）的方法出现使得粒子数目可以进一步增加。到了新千年，CPU主频停止增长，硬件提升的红利也消失了。既然单个核心无法提升，人们转而开始大规模地使用并行的技术。
同时结合了各种加速算法之后，到2005年，模拟粒子的数目已经突破百亿（10^10）。除了对于软件上的改良，在硬件上的加速也被引入数值模拟当中。比如牧野等开发了用于计算引力对的专用芯片Grape。最近，异构的计算平台的使用异军突起，比如在GPU和众核芯片上也取得了非常好的性能和重要的进展。目前在世界最大的一批超级计算机（中国的天河2，美国的泰坦、BG/Q，日本的京）上都运行了超大规模的宇宙学模拟。
2014年之后，最大的数值模拟规模已经使用了超过万亿（10^12）的粒子。如果画出一条模拟数目的历史增长趋势图的话，其速度是超过摩尔定律的。
即便如此高的分辨率还只能够给出宇宙总体的结构，对于我们感兴趣的很多物理现象依然淹没在数值的噪声中了。在无法增加粒子的数目的情况下，局域再模拟（re-simulation）的技术帮助我们大大地拓展了数值分辨率的极限。
通过对一个感兴趣的星系团进行再模拟，分辨率就能够大幅地增加，我们可以清楚的看到原本淹没在数值分辨率下面的子结构在更高分辨能力的模拟中清晰地浮现了出来。很多人可能会怀疑数值模拟的结果是否真的可信。对于一些特殊的系统，数值解的正确性确实有一些分析的方法，但对于高精度的数值模拟来说，类似的分析是很难的。实际上，对相同的初始条件进行多个级别的再模拟，这正是一种数值收敛实验。
也就是低精度的模拟与更高精度的模拟相互符合的部分就是数值上收敛的，通常也是可以相信的部分。有些再模拟最高的分辨率甚至能够达到地球的质量，这对于研究有关暗物质极其相关的性质打开了窗口。
前面我们并没有仔细区分宇宙物质能量组分中的暗物质和重子物质，只是统称为物质。现在我们简单地介绍一下。我们并不知道暗物质的本质是什么。
我们只知道它通常不产生、不吸收、不反射光子，也就是说它不和光发生什么关系，所以叫暗物质。相对于暗物质，我们把更为熟知的恒星、地球、人类、灰尘、空气等等物质统称为重子物质。根据最新的普朗克卫星测量结果，现在人们知道：所有的重子物质占全部宇宙组分的大约4.9%，而暗物质约占到26.8%左右，剩下的是暗能量，约占68.3%左右（误差大概在百分之几）。
既然暗物质几乎不和其他物质发生作用（collisionless），我们是如何知道它的存在的呢？是引力。对于恒星和星系，它们的速度不是很快，会很容易陷在暗物质的势井中无法逃脱，而光的前进路径也会发生偏折，造成所谓的引力透镜。引力势场决定于物质的分布，同时物质的分布是引力系统演化出来的结果。这两者相互依赖，导致不同性质的暗物质，经过去上百亿年的演化之后会产生不同的宇宙结构。
比如运动得更快更“热”的粒子，也就更不容易被束缚，结果产生的结构就会更少。通过对比运行不同暗物质模型的模拟，最终发现如果暗物质过“热”，那将会在定量上背离我们的观测数据。所以今天的标准宇宙学模型中普遍认为存在一个冷暗物质（Cold Dark Matter）的组分。这样的结论很难只从定性的角度分析得到，可以认为这是数值模拟的一大贡献。
基于这一套冷暗物质宇宙学框架，数值模拟能够以非常高的精度给出星系团、重子声波振荡、引力透镜等观测手段的预言，同时可以加入大量从理论上难以计算的效应，估计系统误差，指导未来的观测巡天等等，真正架起了从理论到精确观测之间的桥梁。
前面的介绍主要侧重于N体，N体模拟对于无碰撞的暗物质是非常有效的，但是我们实际观测到的星系，其分布和暗物质的分布之间存在着紧密的联系。打个比方，这就像我们可以通过地球夜晚灯光的分布来了解人口聚集程度一样。星系就是点亮我们宇宙的灯火，我们的数值模拟同样也需要加入星系。
星系形成的物理原理很复杂，其中一个最为重要的因素可能就是气体的冷却和反馈机制之间的平衡。让星系明亮的是其中的恒星。
前面我们说过，宇宙在开始时是一团粥，恒星和星系都是随着宇宙结构的增长而逐渐出现的。恒星是正在进行核燃烧的气体球，其密度是非常高的，宇宙原初的气体从原初极端均匀的状态，塌缩成致密的气体云，进而在其中孕育恒星。但随着气体变得更加致密，引力势井也更深，气体的势能转化成热，这个热能将撑住引力而阻止进一步塌缩称为恒星。但是重子物质（气体）不同于暗物质，它能够冷却，也就意味着可以进一步的塌缩。
实际上在模拟中，通常使用不同类型的粒子来分别代表气体和暗物质。随着恒星的形成，还要加入代表恒星的粒子。
恒星也并不是静止不变的，它们有自己的生命周期。有星风，有超新星爆发。这会把恒星粒子的质量再次交还给气体组分，同时改变了气体的温度、密度、压强等状态，进而改变下一代恒星的形成。也就是所谓“反馈（feedback）机制”。除了恒星，活动星系的中心核区（AGN）也可能提供反馈。这些反馈的影响甚至可以超出星系本身，最终剧烈地改变宇宙中的星系、以及星系中恒星的状态。
最后我们来看更大尺度上的一段视频，来自Illustris Project。左边的盒子显示的暗物质的分布，右边的盒子是同一个模拟的气体成分，盒子的转动是由于视角的转动，而不是宇宙真的在旋转。可以看到右边的气体成分基本跟随着暗物质的引力势井，在塌缩的过程中温度升高，在密集的地方形成星系，45秒之后看到有大量的爆发现象，这些就是气体的反馈机制的效应。
如果宇宙是一个海洋，星系就是一个个岛屿，而宇宙数值模拟就像一个水族馆。在其中我们不能指望发现一个我们没有放入的物种，但是我们可以学习到已知物种的行为的细节和它们之间的相互影响。天文学是无法进行孤立重复实验的，而宇宙数值模拟可以让研究者掌握它全部信息，因此在研究工作中成为了理论和观测之外的重要一极。现在的数值模拟虽然取得了很大的成功，但其中还有很多细节值得商榷。
未来是宇宙学气体数值模拟的大发展期，这不仅反映在计算量上（分辨力的提高），还在于气体演化的物理机制。从第一原理到恒星形成、星系形成之间还存在着一条鸿沟，这既是一个巨大的挑战，也是推动宇宙学和星系形成数值模拟发展的动力。