導論(Introduction)
眾所周知的是二階彈性常數是一個二階四秩張量（cijkl），三階彈性常數是三階六秩張量（Cijklmn）；更高階的彈性常數還包括四階彈性常數（Cijklmnpq），五階彈性常數（Cijklmnpqrs）和六階彈性常數（Cijklmnpqrsuv），分別是八秩，十秩和12秩張量。目前研究報導的最高彈性常數為6階彈性常數。二階彈性常數的計算由於張量分量較少，同時在計算過程中主要是採用了線性hooker定理，所施加的應變數很小，因此在材料力學性能表征中得到了廣泛的套用，三階彈性常數描述了非線性Hooker定理或者非線性力作用下的材料的力學回響問題。三階彈性常數矩陣形式十分複雜，即使對於立方晶體結構也有6個獨立分量，對於對稱性更低的晶體結構則分量更多，如果採用能量－應變的方法來計算各個獨立分量，則計算量相當可觀。如立方晶體（點群Oh，O，Td）有六個獨立分量，則需要六個不同的應變模式，得到六個多元一次方程組，聯立求解得到各個分量數值。若採用應力－應變（Strain－Stress relations）則可以明顯減少應變模式的數量，但主要問題在於這要求第一原理計算軟體具有計算晶體Cauchy應力張量的能力，然而，目前廣泛採用的DFT計算軟體，如VASP，Wine2K等等是不能直接得到Cauchy應力張量的，只能計算特定應變下的應變能數值。Materials Studio軟體中的CASTEP模組是目前為數不多的具有直接計算應變結構Cauchy應力張量的軟體，因此CASTEP模組在計算材料的二階彈性常數方面十分的方便。此外，對於三階彈性常數，目前沒有軟體可以直接計算，需要研究者自行設計方法進行計算。三階彈性常數可以描述材料在高壓下的力學回響情況，鑒於目前高壓物理學，行星結構科學等相關領域的飛速發展，對於超硬材料的非線性力學常數受到了越來越多的關注。採用超聲腔共振法可以方便的測量材料的二階彈性常數，但對於非線性彈性常數，試驗方面進展十分的緩慢，時至今日，大部分超硬材料的高階彈性常數仍然是未知的。
1 廣義Hooker定理以及晶體的Neumann原理
(The generalized Hooker's Law and Neumann Principle)
1.1 廣義Hooker定理
(The Generalized Hooker's Law)
與傳統線彈性力學中廣泛採用的Hooker定理相比，廣義Hooker定理可以認為是包含了高階非線性應力和應變關係項的Hooker定理的Taylor級數展開形式，因此從數學意義上來講這種應力對應變的階數可以無限制的進行，正如前文所說，目前廣泛採用的彈性常數為應力對應變展開的線性項，展開係數即為Cijkl，Cijkl就是彈性常數，由於根據張量運算法則可知Cijkl是一個四秩張量，描述兩個二秩張量應力（stress）和應變（Strain）之間的關係。Cijkl矩陣元素有3^4個，考慮到Lagrange應變為對稱矩陣，同時應變自由能與應變路徑無關，可以用一個6×6的矩陣來描述，有36個矩陣元素。進一步的矩陣元素化簡來自於晶體結構點群對稱性對物理學性質的限制，即Neumann原理。例如對於三斜晶體（Triclinic Crystal Class），獨立矩陣元素為21個，如果是正交晶體（Orthorhombic Crystal Class）則減少到9個，對稱性最高的立方晶系（Cubic Crystal Class），僅為三個（C11， C12和C44）。
Cijkl
Generalized Hooker 's Law簡單歸納一下Cijklmnpq....的性質：
Cijkl: Second order 4th rank tesnor, 6×6 matrix with 36 matrix elements and reduces to 21 independent elements for Triclinic crystal class (The matrix form see Fig);
Second order 4th rank tensor: CijklCijklmn: Third order 6th rank tensor, 21×6 matrix with 126 elements and for triclinic crystal class the independent numer is 56;
Third order 6th rank tensor (Cijklmn)Cijklmnpq: Fourth order 8th rank tensor, 56×6 matrix with 336 elements and for cubic crystal class the independent nuber is 11;
Fourth order 8th rank tensor (Cijklmnpq)Using Voigt notations, we have Cijkl－Cij， Cijklmn－Cijk and Cijklmnpq－Cijkl;
Voigt Notations:
examples: C1111－C11， C111111－C111 and C11111111－C1111；
C1122－C12；C112233－C123； C11112233－C1123；C2332＝C2323＝C44 and etc. etc.
| 11 | 22 | 33 | 23(32) | 13（31） | 12(21) |
| 1 | 2 | 3 | 4 | 5 | 6 |
1.2 Neumann 原理
（Neumann Principle）
Neumann原理指出，任何晶體結構的物理性質所具有的對稱性不低於晶體點群的對稱性，這表明張量分身最少具有晶體點群對稱性，將晶體點群對稱操作作用到各個張量分量Cijkl上，得到新的張量Cmnpq ，則Cijkl＝Cmnpq ，即張量各個分量相等。對於三階和四階彈性常數則可以表示為
Neumann Principle對於高秩張量，上述計算必須在計算機上編程進行，否則計算量是十分巨大的。對於晶體結構，如果張量採用了Voigt標記法之後，張量在晶體點群對稱操作下的變化相當於對張量下標進行變換，對Voigt標記數字進行對稱變換得到的下標分量的對應關係與對張量本身變換得到的分量關係是一致的。現在存在的問題是晶體點群包含的對稱操作的個數很多，如果將每個點群元素對張量進行操作，則不能有效的減少運算量，並且這樣的操作大多數是重複的。根據點群的性質，對於每個點群總能找到幾個代表性的操作元素（Representative Operations），點群中所有的對稱操作元素可以通過代表性元素的冪乘積得到，因此代表性元素稱為點群的生成操作（generators）。因此只要找到點群的生成操作，將生成操作對張量進行變換即可得到點群對稱性限制下獨立矩陣元素。對於32個空間群點群的生成操作在大部分群論方面的書籍中均可以找到，當然也可以藉助群論分析軟體Isobyu獲得。
簡化後的二階，三階以及四階彈性常數矩陣形式依然十分複雜，這裡給出幾個例子：
examples： Orthorhombic crystal class: Space Group Pmm2; Point group mm2
Cijkl （Cij）
Cijkl for Orthorhombic crystalsCijklmn （Cijk）
Third order Cijk for Orthorhomic crystalsCijklmnpq （Cijkl）
The forth order Cijkl for Orthorhombic crystalsCubic Crystal Class: Space Group P-43m; Point group -43m
Cij
Cij for Cubic crystalsCijk
Cijk for Cubic crystalsCijkl
Cijkl for Cubic crystals（下文未註明時彈性常數均為Voigt notation， Cijkl－Cijklmnpq）
2 應變自由能與彈性常數的關係
（The Relationship between Strain energy and Elastic Constants）
應變自由能（Strain energy）定義為應變晶胞與初始平衡結構總能量的差值，應變自由能與彈性常數不能直接對應。實際上如果將應變自由能除以初始晶胞的晶體得到新的參數稱為應變自由能密度 （Strain energy density），該參數對Lagrange應變的Taylor級數展開直接對應彈性常數的某個組合或者彈性常數本身。在應變形式比較複雜的情況下，應變自由能密度展開係數一般為彈性常數的線性組合形式，因此如果要通過應變自由能密度來計算彈性常數則需要求解一組多元一次方程組，方程組的個數即為獨立彈性常數的個數。如對於四階彈性常數Cijkl對於正交晶體有42個獨立變數，則至少應該有42個多元一次方程組，同理立方晶體至少為11個。
Strain energy density and elastic constants3高階彈性常數的獲得與計算
（The calculations for Higher order elastic constants）
3.1 實驗測量
（Experimental determination）
高階彈性常數的試驗測量主要是採用超聲腔共振法（Ultrasonic Resonance Method），如對於二階彈性常數Cij，可以直接採用該方法，設定不同波長的縱波（P Wave）和橫波（S Wave），其中縱波對應晶體結構中的疏密變化，即拉伸模量數值，如C11，C22，C33等；而S波則對應剪下變形即C44，C55，C66等。在施加混合波形情況下可以得到C12，C13等混合應變模式下的彈性常數數值。對於C11而言可以施加一個沿[100]方向傳播的P波得到，具體關係為Cijkl＝pv^2(p: density; v: Velocity)。對於三階彈性常數，可以通過測量二階彈性常數與溫度之間的關係得到，同時注意到Cijk是Cij對應變一階導數，因此也可以通過對晶體預先施加特定的應變（應力模式），然後測定聲速的變化來對曲線進行擬合得到。更高階的彈性常數如Cijkl則需要通過計算Cijk在特定應變下的波速變化計算，由於高階彈性常數獨立變數很多，因此試驗測定高階彈性常數遇到的極大的困難，到目前位置大部分晶體的三階以上彈性常數都是未知的。
3.2 理論計算
（Theoretical Calculations）
由於計算機的飛速發展和固體理論，凝聚態物理理論方面的巨大進度。高階彈性常數目前可以通過理論計算獲得，主要是基於密度泛函理論的第一原理計算（The ab-initio calculations based on density functional theory).正如前文提到的，目前計算方法有兩個：應力－應變關係和應變能－應變關係。前者基於廣義Hooker定理，後者基於應變能密度的Taylor級數展開。兩個方法各有優勢，Stress－Strain方法最大的特點是所需要計算的應變模式很少就能計算出全部的彈性常數，如Cubic晶體，只需要一個應變模式即可得到C11，C12和C44，但若採用應變能密度－應變關係則需要設計三個應變模式分別計算獲得三元一次方程組，求解得到全部常數，應變－應力計算主要需要計算出晶體的Cauchy應力張量，但由於晶體Cauchy應力張量對邊界條件敏感，因此如果精度計算方法對擬合精度都有嚴重影響；應變能密度－應變關係特點是計算量很大，但數值精度穩定，計算結果比較可靠。
4參考文獻
（References）
1. Hiroshi Soma, Yosio Hiki, Journal of the Physical Society of Japan, 37 (2), 544 (1974)
2. Shigeru Mori, Yosio Hiki, Journal of the Physical Society of Japan, 45 (5),1449 (1978)
3. Yosio Hiki, Annual Review of Materials Science, 11, 51-73 (1981)
4. Xanthippi Markenscoff, Journal of Applied Physics, 50 (3), 1325 (1979) (The 4th order elastic constants for hexagonal and trigonal crystal classes are discusssed)
5. P.B. Ghate, Journal of Applied Physics, 35 (2), 337 (1964)
6. T.S.G. Krishnamurty, Acta Crystallgraphy,16, 839 (1963)
7. REMI Brendel, Acta Crystallgraphy, A35, 525-533 (1979) (The independent fourth order elastic constants for all crystal classes are listed)
8. David Y. Chung and Yuan Li, Acta. Cryst, A30, 1 (1974) (The 4th, 5th and 6th order tensors are discussed in this paper and the matrix forms of them are also tabulated, very good article!)
9. Lars Fast and et.al., Physical Review B, 51 (24), 17431 (1995) (The Second order elastic constants for Hexagonal crystals were studied)
10. S.K.R. Patil and et.al., Physical Review B, 73, 104118 (2006) (The second order elastic constants for Orthorhombic and Tetragonal, Cubic crystals)
11. O.H. Nielsen, Physical Review B, 34 (8), 5808 (1986) (both of the Second and third order elastic constants of Diamond were calculated)
12. Jijun Zhao and el.al., Physical Review B, 75, 094105 (2007) (The general method for the calculation of higher order elastic constants using density functional the
