查看原文
其他

费米拓扑已不凡 拓扑玻色若安闲

特约通讯员 知社学术圈 2019-06-30

海归学者发起的公益学术平台

分享信息,整合资源

交流学术,偶尔风月


天 问


离骚一唱两千年

点点滴滴未咏传

渡越时空相望远

方知玻色亦婵娟


1. 引子


这是一篇需要收拾心情、端坐下来阅读的文章。其中的故事是凝聚态物理人工作状态的一种写照:衣带渐宽终不悔、为伊消得人憔悴。


如果您近年来关注过、或正在关注凝聚态物理与材料方面的研究,对“拓扑绝缘体”、“外尔半金属”等名词一定不会陌生,其实您想陌生也不行。这些材料中的电子能带结构,因为其波函数所具有的“拓扑”属性而变得不平凡,进而导致各种奇异的物性,比如绝缘体的“表皮”竟然很导电及具有“手性反常”特性,等等。


那么,什么是能带上的拓扑呢?笔者知道,要将这个问题说到您懂我懂,其实比解决这个问题本身的难度小不了多少。“知社学术”不久前曾经刊发过一篇科普文章,请点击抽丝剥茧是为真一片匠心落子辰。如果您需要在去物理工作室的路上临时买一份早餐,那您御览这篇科普,应该就能饱了。


事实上,三维空间中的拓扑其实不难理解,它体现了物体形状的某些不变性。具体来说,就是允许对物体任意扭曲、变形 (但不允许切割或粘接),物体的某些性质依然不变。这样说来有些抽象,举个例子:一个物体上面“有几个洞”(不是凹进去的那种洞,是前后贯通的那种洞)就是一个拓扑不变的性质,而拓扑不变量就是“洞的数目”。一个带把手的杯子和一个甜甜圈是拓扑等价的,因为它们都只有一个洞。但这个杯子与一个橙子则不等价,因为后者没有洞。由于拓扑不变量是在任意变形下都不变的东西,它天然具有很强的“稳定性”,不会因为形状在细节上的变化而变化。


在变化中寻找不变,这恰恰是物理学家所追求的。只是,在能带体系中,我们考察的对象不再是三维空间中的形状,而是生活在抽象的希尔伯特空间的波函数 (也可以设想倒易空间中的“物体”形状)。在这个无穷维线性空间中,仅仅靠常识的想象力是无法理解其中的拓扑的。能带体系中的拓扑不变量,只能借助代数拓扑学中的工具才能求得。不过,正如前述例子中“洞的数目”一样,所有的拓扑不变量都是整数,而且不因系统的一些细节变化而发生变化,具有高度的稳定性。更妙的是,当一个材料具有非零的拓扑不变量时,它在边界上会呈现出具有特异性的边界态:像是“一方通行”的光波、有着“无限螺旋”样的色散……这些略显中二的物态。图1和视频1 所示即为几个生动的实例,虽然看君欣赏起来也还是要费一番脑子。


图1. 左图:具有两个正的和两个负的Weyl点的拓扑光子晶体之表面态。蓝色红色的圆锥是体态的投影,黄色是表面态的色散。由于螺旋面的性质,在不同能量考察表面态的等能线,会发现它们 (绿色、黑色、蓝色) 随着能量的变化在转动,转动的方向取决于外尔点的正负。来源:Science 359, 1013 (2018)。右图:这是两个狄拉克点投影在表面布里渊区上同一高对称点时,表面态的螺旋面色散。其中黑色圆锥是体态投影,四个彩色的螺旋面分别来自两个狄拉克点。来源:Nature Physics 12, 936 (2016)。


https://v.qq.com/txp/iframe/player.html?vid=k13449ialp6&width=500&height=375&auto=0

https://v.qq.com/txp/iframe/player.html?vid=j13449kuq55&width=500&height=375&auto=0

视频1. 在磁性拓扑光子晶体中模拟出的边界上单向通行的光波。上部是拓扑不变量为0 的情形,可以看到这时候光波可向两边传输。下部是拓扑不变量为1 的情形,此时由于手性被破坏,光在边缘处只能向右传播。注意,每种情况下,晶体上部边缘有个小圆圈,那是偶极天线输入光波的地方。对比两种情况,我们看到,上部视频中,从小圆圈输入的光波可左可右传播;而下部视频中,光波只能向右传播。这里,“跨过障碍”和“一方通行”是一回事。由于没有反向行驶的行波解,所以即便有散射物,也无法把波“弹回去”,因为压根就没有左行波与入射波形成驻波。来源:Nature 461, 772 (2009)。


2. 以 Weyl 之名


请允许我们再稍微物理一点。以外尔半金属为例,里面特别有意思的能带结构是所谓的“外尔点” (Weyl point)。外尔点附近动量空间里的电子结构,可以用一组对每个三维动量都包含两个能级 (能带) 的哈密顿量来描述。将动量坐标的原点取在外尔点位置,以分量的形式展开。在略去各种不重要的系数,也就是将困难和复杂性隐藏之后,系统哈密顿变得简洁、对称、美妙起来:


求解本征值问题后,对每个k,可得到两个能量本征值以及相应的能量本征态。前者就是著名的外尔锥状色散关系 (与狄拉克锥看起来是一样的),如图2a 所示。能量本征态可以用一个含有两个复数的向量来表达,但更直观的做法是画成图:把复数向量表示为“Bloch 球面” (图2b) 上的一个点 (图2b 中的箭头)。其中,球的“北极”是 (1  0) 这个态,“南极”是 (0  1) 这个态,“纬度”  θ 表示向量的南北极分量之间的“权重”比例,而“经度”φ 则表示南北极分量之间的相对 (复数) 相位。


图2. (a):外尔锥形式的色散。横轴表示三维动量空间中沿任意方向到外尔点的距离。虚线与外尔锥的交点对应于 (c) 中的等能量面。(b):用于表示本征态的Bloch 球面,摘自Wikipedia。(c):等能量面上的本征态,在图中用末端位于Bloch 球面上的红色箭头表示。


经过简单 (其实也不简单) 计算,不难发现计算结果中有一个特点:随手在动量空间中画一个包裹住外尔点的闭合曲面 (不妨将其取为一个球面,比如图2c 中所画的外尔锥上一个等能量面),如果把面上各个动量点处能量较低的本征态在Bloch 球面上的点都描出来,就不多不少,刚好铺满了整个Bloch 球面。这就是物理的精妙!


奇迹由此产生——这样的外尔点,是不能凭空被“消灭”的。如果我们随手画的闭合曲面并不包含外尔点,同样执行后面的步骤后,就会发现:这时在Bloch 球面上描出的点,只会盖住球面的一部分。极端地讲,哪怕允许这个被盖住的Bloch 球面部分被任意拉伸,最多也只能“几乎”铺满整个球面 —— 会留下一个洞。而“拓扑”,正如第一小节所言语,笼统地就是关于几何图形中“洞”的数学。


正是上述波函数结构“有”或者“没有”洞的差别,使得我们无法凭空把外尔点式的能带相交解除。因此,我们说,外尔点是受拓扑保护的能带交点。


我们看到,当一个动量空间球面上的波函数具有非平凡的拓扑结构时,这个球面一定会包裹一个拓扑能带交点。这就是外尔点或者说是“拓扑半金属”的成因。而当整个动量空间某个能带的波函数作为整体具有非平凡的拓扑结构时,这个能带与其它能带之间的空白,即所谓的“能隙”,就具有了非平凡的拓扑性质。这就是拓扑绝缘体的来源。拓扑半金属和拓扑绝缘体是拓扑能带理论中的两类最重要的体系,也是研究结果最丰富的体系。


能带上的拓扑之所以让人着迷,不仅因为它能够带来“解不开”的能带相交。波函数的全局结构,还可以导致电子在材料体内具有包括几何相位在内的非平庸特性,以及在材料表面 (或边界) 上的奇异状态。比如,一个重要的特征是外尔半金属表面上电子的“费米弧”。


从这个简单例子,我们还可以感受到一点拓扑能带理论的精髓:矩阵形式的代数给出了从动量空间到抽象向量空间 (希尔伯特空间) 的一个映射,而拓扑则描述了映射关系中的“洞” (或者没有“洞”)。到此,看君应该感受到了:物理人为了将他们心中的“物理”表述给您,是绕了万水千山的!^_^



3. 玻色拓扑


好奇的人会进一步问:由矩阵描述的物理系统有很多,而只要有空间上的周期性 (晶体),动量空间的概念也能广泛适用于其它粒子或者波动。那么,我们可不可以超越电子来谈这些概念?对其它粒子,是不是也可以如法炮制这些基本物理?


没错!在前两节的论述中,我们并没有强加微观粒子必须是电子或费米子这一要求。一个基本粒子,要么是费米子,要么是玻色子。事实上,玻色子体系也可以具有能带上的拓扑结构,甚至经典的波动也可以。比如在所谓的光子晶体和声子晶体这些人工结构中,拓扑非平庸的能带结构都已经被发现了 [1, 2]。


然而,如何在天然晶体材料中找到拓扑非平庸的、真正意义上的玻色子能带结构,却是一个比较难的问题。


首先,选择目标:天然晶体中,可以看着是玻色子的有两类效应:晶格振动、自旋激发。它们是晶体中常见的玻色子 (准粒子)。


其次,工欲善其事、必先利其器!为了攻克这一问题,要选好用来做探测实验的工具。如何测量晶体材料中的玻色子能带?非弹性的中子散射似乎提供了理想的工具,因为无论是晶格振动、还是自旋激发,都可以被中子散射直接探测到。


果然,世上英雄无数,使得这个时代新鲜事物已经不多。您我所能想到的事,很多早已经发生过了。2015 年和2017 年,由美国MIT 和英国ISIS 等单位组成的团队,就分别使用非弹性中子散射对玻色子体系中的“拓扑能隙”进行了探索。其中所涉及到的玻色子是携带有自旋量子数1 的自旋波激发 [3, 4]。通过比对散射实验和理论计算,他们认为实验结果很好地符合了理论预言的谱学形貌 (图3)。而理论告诉我们,这个能带结构是拓扑的。所以,实验也就似乎圆满地完成了任务 —— 玻色子拓扑绝缘体已经找到。


图3. Cu(1,3-bdc)在7 特斯拉磁场中的粉末样品非弹性中子散射数据 (左图) 和相应的理论计算 (右图)。根据计算,图中能量最高和最低的能带具有非零陈数 (分别为 +1 和 -1 ),因而是拓扑非平庸的。图片摘自文献 [3]。


4. 桃花流水


细心的人也许会产生一个疑问:实验中到底看到了什么,就让 MIT 和 ISIS 这帮人“轻而易举”地宣称达到了预期目标?


图3 中,我们主要看到的是一系列正常的能隙和拓扑绝缘体中同样具有“反带”性质的能隙。换句话说,最主要的实验证据就是观测到的能谱形状看起来和理论预言的拓扑能隙附近的能谱形状很相似。可是如果没有理论计算的帮助,如何知道能隙的两侧是“反带”呢?


答案有些让人失望:其实还真是没法知道!!!


这里有一些细微亦博大之处。首先要问:人们对拓扑绝缘体的理解到底是否正确,拓扑绝缘体当年又是如何被实验证实的呢?我们知道,体能带中的非平庸拓扑,通过所谓的“体-表全息关系”,会造成独特的表面态。因此,如果实验上观察到受拓扑保护的表面态,那当然可以为拓扑能带理论提供强有力的支持。拓扑绝缘体是幸运的,因为用来测量电子能带结构的角分辨光电子谱 (ARPES) 技术,恰好也对表面电子结构高度敏感,因此拓扑绝缘体的表面态 (具有狄拉克锥形式的色散) 很快在实验中被观测到 [5]。一切都很好!


可是,看君注意了:我们是谈玻色子,不是费米子。此时,ARPES 技术就有点捉襟见肘,因为它探测玻色子并不那么来事。遗憾的是,此时中子散射的短板也出现了:中子可以很容易穿透厚重的水泥墙,所以一般用来做非弹性中子散射的实验样品,最少也要有好几克的质量。如果您的样品只是一个“二维”表面,那中子穿过去后带走的信息微乎其微。因此,对仅仅在原子级厚度的“表皮”内存在的表面态,那一点样品量根本不够测,除非您将成千上万的这种“表皮”叠加在一起。可是您要是揭掉“表皮”,它马上就不再是有表面态的“表皮”了。这就是“皮之不存毛将焉附”的意思,您懂的!


所以我们要问:脱离了对拓扑表面态的观测,是否还能说美国MIT 和英国ISIS 团队所做的中子散射实验就真的找到了玻色子拓扑绝缘体?哦,这恐怕只能说是一个见仁见智的问题。事实上,它就是一个大大的问题。


5. 山重水复


最近,事情又有了新的突破。既然用中子散射看“表皮”没戏,而我们又不能不用中子散射,那就得放弃“表皮”,转而深入到块体内部去找机会。由此,物理人的思路就从“拓扑能隙”转向了寻找玻色子中的“拓扑能带交点” (或者说“拓扑半金属”态)。这是关键点,这是在说体系内部的事、三维大块体系的事、体能带的事。您看,机会来了。


注意,拓扑绝缘体的体能带看起来就和普通绝缘体一样 (有能隙),需要借助理论计算或对拓扑表面态的观测,才能判定能带的拓扑性。与此对比,拓扑半金属的体态具有拓扑保护的能带相交,例如图2 中的外尔点。当然,我们说仅仅靠能谱形状是无法确定“拓扑能隙”的,但与此不同,理论上我们知道:如果从一个能带交点附近沿着所有方向看能带的色散,如果所有色散都是线性的,那该交点必然是拓扑非平凡的。这一特征具有良好的可观测性,对中子散射来说无疑是个好消息:体能带中有了与拓扑相关的可观测特征。


2016 年,复旦大学陈钢团队 [6] 和德国马普学会的Alexander Mook 等 [7] 率先在某些非中心对称及具有铁磁性的晶体中预言了磁振子外尔点的存在。不过,满足这种条件的材料很少,这些理论至今还没有获得直接的实验验证。


一年以后,中科院物理所方辰研究组和北京大学量子材料科学中心李源研究组在计算合作中发现,磁振子的拓扑能带相交这一特征竟然存在于一大类常见的反铁磁体系中 [8]。之所以说“常见”,是因为这一理论预言对体系对称性的要求很低:只要具有PT (中心反转和时间反演的联合操作) 对称性就够了,而大部分已知的反铁磁材料都满足这一要求。


为了分析方便,他们先是假设体系同时还具有全局的自旋旋转 U(1) 对称性,基于此证明了磁振子能带原则上可以有重叠在一起、却不会彼此湮灭的外尔点,也就是狄拉克点 (图4a)。然后,他们再把 U(1) 对称性去掉,进一步证明此时狄拉克点会张开成为一段 (或者一圈) 节点线 (图4b)。这些拓扑能带相交的方式,在电子体系中还不曾有过具体到材料的预言或实验观测。


图4. (a):对Cu3TeO使用 J-J(J2 =0.134 J> 0) 模型计算得到的磁振子能带结构。D1 - D3 表示狄拉克点。(b):当自旋旋转的 U(1) 对称性被破坏时,位于布里渊区 P 点处的狄拉克点张开成为节点线。(c):Cu3TeO的反铁磁基态和由最近邻相互作用 J构成的网络。摘自文献 [8, 11]。


具体的计算是针对一种化学式为Cu3TeO的反铁磁材料进行的。这个立方晶系体系此前被研究人员起了一个独特的名字:“自旋网” (spin-web) [9]。原因是,如果仅考虑携带自旋 1/2的 Cu 离子之间最近邻相互作用,会发现每个Cu 只有4 个邻居。也就是说,自旋格子虽然是三维的,但其中相互作用“连线”的密度却与二维方格子类似,就像一块稀疏的网筛。事实上,最初李源和方辰研究组之所以会注意到 Cu3TeO6,正是因为其稀疏的自旋网结构 —— 他们原本打算通过掺杂等方式,破坏掉其中貌似“脆弱”的反铁磁序,进而得到各种可能的演生现象,比如高温超导。但事与愿违,随后研究表明这个材料里的反铁磁序一点都不脆弱。但幸运的是,这对能带拓扑的研究却是一个好消息。反铁磁的基态越是稳定,磁振子作为准粒子的寿命就越长,从而也越能给出清晰的谱学信号。此外,Cu3TeO6  结构的高对称性及含有多达12 个Cu 原子的磁原胞,都给实验研究带来了不小的帮助。


6. 春风渡越


在理论预言的鼓舞下,实验的研究得以紧锣密鼓地展开。追逐这一目标的不仅有李源和方辰自己的研究组,由南京大学温锦生、于顺利、万贤纲 (电子能带中外尔费米子的主要发现人之一)、李建新带领的研究团队也迅速投入了工作。按照方辰的说法,“能够被两个实验分别进行检验的理论,是幸运的”。


两个研究团队采取的思路不约而同地相似——使用基于飞行时间原理的非弹性中子散射谱仪,对磁振子能带结构进行完整的四维 (三维动量 + 能量) 探测。这些实验分别在日本的 J-PARC 中子源和美国的橡树岭国家实验室完成。该方法最大限度地运用了中子散射的优势:实验过程中收集到的大数据 (多达几十Gb ) 包含许多个布里渊区中的散射信号,不仅可以用来分析磁振子在布里渊区上带有周期性的色散关系,信号强度本身也直接与磁振子的本征波函数有关。换言之,如果理论的建模足够准确,通过计算给出的本征波函数信息也将能够得到实验检验。


需要强调的是,直接对体态波函数进行测定,这在以往的拓扑能带结构研究中并没有先例 (ARPES 实验可以一定程度上测定表面态的电子波函数,比如电子的自旋 - 动量锁定关系 [10])。


图5. (a):Cu3TeO的第一布里渊区和高对称点。(b) 与(c):中子散射实验得到的沿动量空间[001] 和 [111] 方向的磁激发谱。其中的白线为理论计算给出的色散关系。内插图表示各自图中黑色虚线框内交点附近的色散。摘自文献 [11]。


现在,让我们通过两张图里的实验结果,来看看拓扑能带理论在这个磁振子体系中被检验到什么程度。图5 是来自南京大学团队的工作。在布里渊区的高对称点位置,例如图5b 的 H 点以及 Γ 点,能带的交点清晰可见。对称性分析表明,这些交点具有稳定的拓扑属性。基于线性自旋波理论,研究者们采用一个以最近邻反铁磁相互作用 (图3c) 为主的模型进行计算,较好地重现出实验观测到的散射信号位置。理论计算得到的能带色散关系由图5b、5c 中的白线所示。分析表明,在图5c 中,布里渊区 P 点上的不同能量位置存在三个狄拉克点,靠近这些狄拉克点的磁振子激发可以用狄拉克方程描述。因此,这些准粒子被称为狄拉克磁振子。该结果跟早前的理论预言 [8] 吻合得很好。

 

此外,实验和计算结果表明在图5b 和5c 的 H 和H′位置上,都分别在两个能量上出现了磁振子的三重简并。理论计算表明,在每一个三重简并点的附近,磁振子能带由两条线性色散能带和一条平带交叉组成。因此,这些磁振子为三分量磁振子。这些高对称点上的线性能带交点,包括狄拉克点和三重简并点,均不依赖于理论模型的参数细节,受到材料本身的对称性保护。


图6. (a) 与(b):以 Γ = (1, 1, 2) 为中心的布里渊区内各条高对称线上的磁振子散射信号和相应的理论计算 (从蓝到红信号逐渐增强)。黑色细线表示理论计算得到的色散关系。图(a) 中的橙色和紫色虚线框分别标出了图 (c) 和图 (d-e) 中的高分辨数据的采集区域。(c):P 点能量最高的两支 (各为二重简并) 能带之间的狄拉克点附近的散射信号。数据采集自 P(3/2,1/2, 3/2) 点和附近多个N 点所在的动量平面,并经过以 P 为中心的对称化处理。(d) 与(e):一条H – P – H 路径上的散射信号强度分布和相应的理论计算。黑色细线表示理论计算得到的色散关系,紫色虚线表示最高能量狄拉克点附近的关于P点不对称的信号强度包络。


图6 来自北京大学和中科院物理所团队的工作。实验数据同样分别在P 点和H 点显示出了磁振子的狄拉克点和三重简并点的存在。为了进一步做到对本征波函数的检验,研究者们对体系中的自旋相互作用进行了几近吹毛求疵的建模努力。检验最终的模型是否有效的依据,是中子散射信号强度在四维空间里的细节分布。图6a 和6b 给出了其中一个布里渊区中实验和最终计算结果的比较,两者的吻合度非常的高。


为了做到这一点,该团队使用了一个与南京大学团队略有不同的模型。该模型中,除了最近邻相互作用以外,“第九近邻”的Cu 之间反铁磁相互作用也十分可观,以至于最终的模型里每个Cu 其实具有多达8 个相互作用较强的伙伴。它们的存在也从另一个角度解释了体系中反铁磁序的高度稳定性。在这个高度精确的模型指引下,研究者们对P 点位置上的其中一个狄拉克点进行了“近距离”的高分辨探测,最终在由P 点和N 点构成的动量平面、能量跨度仅有1.2 meV 的狭小范围里,直接观测到了狄拉克锥式的色散关系(图6c)。不仅如此,他们还在H – P – H 这条路径上细致研究了散射信号在狄拉克点周围的分布。尽管路径上的色散关系注定会关于P 点对称,计算表明散射信号将呈现出一个关于P 点不对称的喇叭形包络,如图6e 中的紫色虚线所示,而这也完全得到了实验数据的支持(图6d)。


7. 春华秋实


至此,可以放心地说,基于自旋波理论的预言已经得到了中子散射实验全面验证。这两项实验工作已于本月早些时候分别在Nature CommunicationsNature Physics 在线发表 [11, 12],文章链接也在参考文献区域给出。实验直接验证了拓扑能带理论在玻色子体系中的有效性,发现了一种新型狄拉克点,以及狄拉克和三重简并磁振子共存的新颖拓扑态,加深了人们对于拓扑能带理论的理解。


毫无疑问,玻色子的拓扑半金属态终于在一个天然的三维晶体材料中被直接观测到。这样一项工作不但丰富了材料科学,也为运用中子散射的技术优势推进拓扑量子领域的发展树立了一个非常成功的案例。


此外,尽管两项实验工作都没有涉及磁振子拓扑表面态和拓扑体态输运性质的测量,但两篇文章在讨论章节里都对这些未来的发展方向作出了展望。维度与对称性,是拓扑能带理论的两大立足点,而描述磁性材料对称性的磁群 (1651个) 又比晶体空间群 (230个) 丰富得多。在磁性材料里,以磁振子为代表的玻色子能带拓扑,究竟还会带来怎样的惊喜与洞见,我们拭目以待。


这一成果很自然会带动更多的研究者进行新的探索。而笔者认为,取得这一成果的过程也多方位展现了物理人开展前沿研究的思路、立足点、方法论与进取姿态。我们说“渡越时空相望远、方知玻色亦婵娟”,不只说婵娟是美貌、月光、草木,更是说物理人渡越时空的景象。


最后,作为对纳税人的致敬,我们愿意指出:本文提到的由南京大学、北京大学、中科院物理所研究团队完成的研究工作得到了来自国家自然科学基金委、科技部、中央高校基本科研业务费专项资金支持,也得到中国科学院战略先导项目等多个项目经费的资助。



参考文献

1. Ling Lu, John D. Joannopoulos, and Marin Soljačić, Nature Photonics 8, 821 (2014).

2. Sebastian D. Huber, Nature Phys. 12, 621 (2016).

3. R. Chisnell et al., Phys. Rev. Lett. 115, 147201 (2015).

4. P. A. McClarty et al., Nature Phys. 13, 736 (2017).

5. M. Z. Hasan and C. L. Kane, Rev. Mod. Phys. 82,3045 (2010).

6. F. -Y. Li et al., Nature Commun. 7, 12691 (2016).

7. Alexander Mook, Jürgen Henk, and Ingrid Mertig, Phys. Rev. Lett. 117, 157204 (2016).

8. K. Li et al., Phys. Rev. Lett. 119, 247202 (2017).

9. K. Y. Choi et al., J. Phys. Condens. Matter 20, 505214 (2008).

10. D. Hsieh et al., Science 323, 919 (2009).

11. S. Bao et al., Nature Commun. 9, 2591 (2018). https://www.nature.com/articles/s41467-018-05054-2  

12. W. Yao et al., Nature Phys. (2018).  https://www.nature.com/articles/s41567-018-0213-x 


备注:

(1)    封面描绘的是量子拓扑能带结构的一幅神来之笔,来自https://phys.org/news/

(2)    题头小诗为Ising 所加

扩展阅读

 

铁磁半导体 花落两维里

PRL : 位错的极化和挠曲电效应

中国团队发现镜像狄拉克锥

超导维度、超越维度

中国热电团队三维电子/二维声子实现高ZT

本文系网易新闻·网易号“各有态度”特色内容

媒体转载联系授权请看下方

    您可能也对以下帖子感兴趣

    文章有问题?点此查看未经处理的缓存