Comments
Transcript
複雑構造材料の特性解析グループ - Department of Precision
複雑構造材料の特性解析グループ 京都大学‐ミシガン大学‐フライブルク大学 MicRO アライアンス 工学研究科マイクロエンジニアリング専攻 田畑 修 Abstract: The objective of this alliance between three Universities, Kyoto University, The University of Michigan and Freiburg University, so called MicRO (Micigan, fReibuRg, kyOtO) alliance named after each Universities name, is the establishment of new design and fabrication paradigms for micro electro mechanical systems (MEMS) of next-generation in which nanostructures are integrated with MEMS. In 2005, the 2nd symposium was held on 14th October at The University of Michigan with about 40 participants. Also as a great step for the close cooperation, a student exchange agreement which including tuition waiver was concluded between Institute of Micro Technology, Freiburg University and Graduate School of Engineering, Kyoto University. Key words: MEMS, Michigan, Freiburg, alliance, education 1. はじめに 2004年11月,京都大学,ミシガン大学,フライブルク大学の3大学は次世代微小電気機械融合システム (MEMS: Micro Electro Mechanical Systems)分野の基盤技術の系統的構築を目指して,3大学の学術交 流協定を締結した.京都大学では多くの大学と学術交流協定が締結されているが,3大学による学術交流 協定の締結は初めてのケースである.この3大学はそれぞれアジア,欧州,北米におけるMEMS分野の研究 拠点として活躍しており,それぞれの大学の名称とこの研究分野からマイクロ(MicRO:Micigan, fReibuRg, kyOtO)アライアンスと名付けられた.現在この学術交換協定に基づき,21世紀COE「動的機 能機械システムの数理モデルと設計論」による支援を受けて,毎年シンポジウムが開催されている.本 稿では平成17年度のMicROアライアンス活動内容について述べる.MicROアライアンスの目的と目標の詳 細,およびミシガン大学,フライブルク大学の概要と特徴は昨年度の報告書を参照して頂きたい. 2. シンポジウム 2005 年 11 月に京都大学で開催された第 1 回シンポジウムに引き続いて,第 2 回シンポジウムが約 40 名の参加者を得て 2006 年 10 月 14 日ミシガン大学アナーバーキャンパス,Lurie Engineering Center (写真 1)の GM Conference Room(写真 2)において開催された.京都大学からは教員 6 名,学生 5 名 の計 11 名,フライブルク大学からは教員 2 名,学生 6 名の計 8 名が参加した.シンポジウム前日の晩に は歓迎のレセプションが開催された(写真 3) .シンポジウム当日はミシガン大学工学研究科長 Prof. Gibala(写真 4)の挨拶に引き続いて集積化ワイヤレスマイクロシステム研究センター(WIMS:Wireless Integrated Micro Systems)のセンター長 Prof. Wise,京都大学大学院工学研究科マイクロエンジニア リング専攻の小寺教授,フライブルク大学マイクロシステム工学科の Prof. Zengerle による基調講演が 行われた. 引き続いて 3 大学の教員による MEMS 分野での研究と教育に関する成果報告と 3 大学の学生に よるポスター発表 (写真 5) が行われ, 活発な討論が繰り広げられた. 午後からは Solid State Electronics Lab.などの研究室見学(写真 6)が行われた.また運営委員会が開催され今後のアライアンス活動につ いての意見交換を行った.京大からの参加学生には,シンポジウム開催期間前後にミシガン大学に滞在 してミシガン大学学生との交流を深めることを強く奨励した.その結果,2 名の学生が開催数日前から ミシガン大学に滞在し,研究室訪問や授業聴講を行うという貴重な経験をした.第 3 回シンポジウムは 2006 年 7 月 7 日,フライブルク大学において開催される予定である. 3. 学生交換交流協定 平成 18 年 2 月に,京都大学とフライブルク大学との間で学費不徴収に関する取り決めを含む学生交換交流 協定が締結された.毎年それぞれの大学から 2 名程度の学生が協定大学において学費を支払うことなく学ぶ ことが可能となった.早速この協定を利用してフライブルク大学から学生(Daniel Schifferdecker 君)が京都大学 に半年滞在する.ミシガン大学の学生からも京大での 1 年程度の滞在を希望する旨のコンタクトがあり,現在具 体的な研究テーマなどについて調整中である.ミシガン大学との学費不徴収に関する取り決めを含む学生交 換交流協定の締結は双方の協定案に隔たりが大きいため調整が難航しているが,実現に向けた努力を粘り強 く続けていく予定である.フライブルク大学とミシガン大学との間でも学生の交流が始まっており,MicRO アライ アンスの成果が着実に現れている.残念な事は,京大から協定大学へ行こうとする学生がまだ現れないことで ある.京大の学生諸君もこの MicRO アライアンス協定を利用して積極的に研究の舞台を世界に広げてもらい たいと熱望している. 写真1 Lurie Engineering Center 写真3 レセプション 写真5 ポスタセッション 写真2 GM Conference Room 写真4 Prof. Gibalaミシガン大学工学研究科長 写真6 研究室見学 骨の再生/機能的適応の複雑系モデルと マルチスケールバイオメカニクス 工学研究科機械理工学専攻 安達 泰治 Abstract: As the mechanosensor cells, osteocytes are embedded in bone matrix, making a three-dimentional network system, and are considered to play important roles in the regulatory mechanism of adaptive bone remodeling. They are connected each other via cell processes in the canaliculi, and are believed to sense the mechanical stimulus due to the interstitial fluid flow in the lacuno-canalicular system. This study aimed to clarify the relationship between morphological characteristics of the lacuno-canalicular system and the mechanical condition at the cellular process level in the trabecular remodeling. First, three-dimensional structure of the lacuno-canalicular system in trabeculae was fluorescently visualized and observed using a confocal laser scanning microscope. Relationship between the canalicular orientation determined by fitting the fabric ellipse and the normal vector field derived from the distance from trabecular surface were quantitatively evaluated. As a result, it was found that the canalicular orientation was approximately coincided with the direction of normal vector field that corresponded to the vertical direction to the trabecular surface. Second, a rate equation of trabecular remodeling was proposed considering the morphological characteristics of the osoteocyte network system. Remodeling simulation conducted for a single trabecula under uniaxial compressive loading revealed that the anisotropy of the system was one of the important factors that affects the trabecular remodeling. Key words: Bone remodeling, Functional adaptation, Osteocyte network system, Multi-scale biomechanics 1. はじめに 骨基質中に存在する骨細胞は,力学刺激の感知センサーとして働き,その情報を周囲の細胞に伝達するこ とで,リモデリング活動を調整していると考えられている.骨基質内部に三次元的に分布する骨細胞は,骨細 管を通して細胞突起を互いに結合し,複雑なネットワーク構造を形成している(1).骨のリモデリングによる機能 的適応の過程において,この骨細胞は,基質に負荷される力学刺激を感知し,その情報を細胞間ネットワーク を介して相互に伝達することで,骨形成・骨吸収活動の調整に寄与していると考えられている(2).骨細胞は,こ のようなメカノセンサーとしての機能を発揮する際,骨に生じる応力やひずみの勾配により生じた骨細管-骨 小腔系内の間質液の流れを感知しているという考え方が提案されている(3, 4). 著者ら(5)は,骨基質に荷重が作用した際に生じる骨細管内の間質液の流れ,骨細胞によるせん断応力の感 知,刺激情報の周囲への伝達等から,破骨細胞・骨芽細胞による骨吸収・形成に至る一連の過程を表現する 骨梁リモデリングの数理モデルを提案した.また,骨細胞が骨基質内において骨細管中に伸ばす細胞突起の 配向が異方性を有することから,骨基質の変形に伴う間質液の流れは,骨細管の配向性や形状に強く依存す る可能性を指摘した(6). 本研究では,骨細胞が形成する細胞間ネットワークに関連したメカノトランスダクション過程に着目し,骨小 腔-骨細管系の構造と細胞レベルでの力学状態との関連性を明らかにすることを目的とした.まず,実際の海 綿骨の骨梁内部における骨小腔-骨細管ネットワーク構造の観察を行い,その形態的特徴について検討する. 次に,これまで構築した数理モデル(5-8)において,特に骨細管の配向性に着目し,単純な骨梁モデルを対象と したリモデリングシミュレーションにより,骨細管-骨小腔系構造形態の配向性が骨リモデリング現象にもたら す影響について考察する. 2. 骨細胞ネットワークの観察手法 2.1 骨細管・骨小腔の蛍光観察 試料として,豚(生後 8 ヵ月)の脛骨近位部より,Fig.1 に 示す 3 mm×3 mm の角柱状の骨組織片を切り出した.ま ず,骨髄を除去し,70%エタノールに 2 日間浸漬し,さら に,暗所にて 50 μM の RH414(膜電位感受性色素)に 4 日間浸漬した.その後,角柱状骨組織片の表面に位置す る骨梁の中から,脛骨長軸方向(荷重方向:Fig.1 の矢印 方向)から見て,円形となる骨梁断面(r-θ 断面),および, それに垂直な骨梁断面(z-r 断面)の 2 箇所を任意に選び 出した.蛍光観察には共焦点レーザー走査型蛍光顕微 鏡(Zeiss, LSM-510)を用い,骨細管および骨小腔の蛍光 画像を 8 bit の輝度データとして PC に記録した. Load direction Cancellous bone Cortical bone Fig. 1: Cancellous bone specimen excised from swine proximal tibia 2.2 蛍光観察画像から骨小腔・骨細管形状の抽出手法 前節の手順により得られる蛍光観察画像は,骨細管内部と骨小腔内部が共に高い輝度を持つ.骨小腔の直 径は,数 10 μm 程度であるのに対して,骨細管の直径は,数 100 nm 程度と小さい.そこで,本研究では,Fig.2 に示す以下の手法により,蛍光観察画像の輝度データから,骨細管および骨小腔それぞれの 3 次元形状を抽 出する.まず,骨梁断面の蛍光観察画像から得られた任意の座標 x における輝度値を L(x)とし,閾値 Lt を用い て,L(x)> Lt で定義される閉領域 D1 を求める.このとき,D1 には,骨小腔部と骨細管の両方が含まれている (Fig.2(a)).次に,細胞突起の半径よりも大きな値δ( >0)を用いて,D1 表面から距離δ だけ内側の領域 D2 を導 出する(Fig.2(b)).さらに,D2 から距離δ 外側の領域 D3 を導出する(Fig.2(c)).これら一連の作業により導出さ れた領域 D3 は,骨細管のような細い突出部を含まない領域として定義される.本研究では,この領域 D3 を骨 小腔部と定義し,領域 D1 が領域 D3 により切り取られる領域 D4 を骨細管部とする(Fig.2(d)).なお,領域 D2,領 域 D3 を導出する際,境界面の形態を厚みδ 分変化させる手法として,Level Set 法(9)を用いる. −δ +δ (a) Domain D1: L(x) > Lt (c) Domain D3 = D2 + δ (d) Domain D4= D1 − D3 (b) Domain D2 = D1 −δ (Lacunae and canalicula) (Canaliculi) (Lacuna) Fig.2: Extraction of lacunae and canaliculi from fluorescent image 2.3 法線ベクトル場 n(x) 骨細管の配向性と骨梁の形態的特徴とを関連付けるための指標として,法線ベクトル場 n(x)を用いる.まず, 対象領域内の任意の点 x から骨梁表面までの最短距離を d0(x)としたとき,符号付距離関数 d(x)を骨梁内部内 部を 負, 骨梁表面を 0 , お よ び 骨梁外部を 正と し て 定義す る . さ ら に , d(x) の 単位勾配ベ ク ト ル n( x ) = ∇d ( x ) / | ∇d ( x ) | により,法線ベクトル場n(x)を定義する.このようにして定義された n(x)は,骨梁表面か らの符号付距離関数 d(x)の等値面に対して,常に垂直な単位ベクトルとなる. 3. 骨細胞ネットワークの観察結果と考察 骨梁の蛍光観察画像の輝度値から,骨梁内の骨細管-骨小腔ネットワークの3次元イメージベーストモデル を構築した.その r-θ 断面,および z-r 断面を Fig.3(a),(b)に示す.なお,Fig.3(a)における右上隅部,右下隅部, 左下隅部,ならびに Fig.3(b)における右上隅部において高い輝度値を持つ領域が見られるが,これらの領域 は,骨梁の外部である.Fig.3 に示すように,骨梁の r-θ 断面では,骨細管は,骨梁中心から放射状に伸展して おり,また,z-r 断面では,骨梁表面に対して,垂直な方向に平行に並んでいる様子が観察された. 次に,第 2.2 節の手順により抽出された骨細管形状および第 2.3 節において定義した法線ベクトル場 n(x)を Fig.4, および Fig.5 にそれぞれ示す.また,Fig.4 の各断面において,骨梁表面から約 40 μm 内側の点を任意 に選び,それぞれ xr-θ,xz-r とし,それらの点を中心とする 15.5 μm×15.5 μm の正方形領域をそれぞれ Region A,Region B とする.さらに,xr-θ,x z-r の各点における法線ベクトル n(xr-θ),n(xz-r),および,各領域内の骨細管を ファブリック楕円体を用いて近似的に表現したものを Fig.6 に示す.Fig.4 と Fig.5 を比較すると,r-θ および z-rの 両断面で,骨細管の配向方向と法線ベクトル n(x)がほぼ一致していることがわかる.さらに,Fig.6 において,法 線ベクトル n(x)とファブリック楕円体の長軸とのなす角をθNF としたとき,Fig.6(a),(b)におけるθNF は,それぞれ, 1.51, 4.83 deg と小さな値を示した.ファブリック楕円の長軸方向は,骨細管の配向方向を示すことから,骨細管 の配向方向は,法線ベクトル n(x)に一致する傾向があると考えられる. 6.0 4.0 2.0 0.0 Region A -2.0 -4.0 Fabric ellipse Normal vector direction -6.0 -6.0 -4.0 -2.0 0.0 (a) r-θ section (a) r-θ section (a) r-θ section 2.0 4.0 6.0 (a) Region A 6.0 4.0 2.0 0.0 -2.0 Region B -4.0 Fabric ellipse Normal vector direction -6.0 -6.0 -4.0 -2.0 0.0 (b) z-r section Fig. 3: 3D reconstruction of lacuno-canalicular system in trabeculae (b) z-r section Fig. 4: Extraction of canaliculi morphology from fluorescent image (b) z-r section Fig. 5: Normal vector field n(r) and the signed distance function d(r) [μm] 2.0 4.0 6.0 (b) Region B Fig. 6: Normal vector and fabric ellipse of canaliculi [μm] 4. メカノセンサーネットワークの数理モデル 4.1 骨細胞ネットワークを考慮した骨形成・骨吸収則 骨細胞間ネットワークによるメカノセンシングを考慮した骨リモデリングメカニズムの一つとして提案したモデ ルの概略は,次の通りである.まず,骨に力学的負荷が作用すると,骨基質に変形が生じ,骨細管および骨小 腔の内部を満たす間質液に圧力勾配が生じる.これにより,骨細管内で間質液が流れ,骨細胞の突起表面に せん断応力が作用する.骨細胞は,これを力学刺激として感知し,その情報を,細胞間ネットワークを介して骨 表面の骨形成・骨吸収の役割を担う骨芽細胞や破骨細胞へと伝達する.前報(10)では,均質化された骨組織に マクロスケールに与えられた力学的境界条件から,ミクロスケールにおいて間質液流れにより細胞突起上に作 用するせん断応力を求め,さらに再びマクロスケールに戻って骨細胞による刺激感知・情報伝達から骨梁表面 移動を考える骨リモデリング現象の数理モデルの枠組みを構築した.本報では,この数理モデルに対して,骨 細管の配向を考慮し,その特性を抽出するため,数理モデルを単純化した. 4.2 数理モデルの単純化 骨細胞が感知する力学刺激量を定義する.骨梁内の任意の点 xb において,単位体積中に N 本細胞突起が 存在するとき,i 番目の細胞突起の伸長方向単位ベクトル npi に対して,骨細胞が感知する力学刺激 Soc を S oc ( x b ) = ∑ | ∇P ⋅ n pi | (3) N と定義した.この力学刺激を細胞間ネットワークを介して骨表面に存在する細胞に伝達する際,細胞間のコミュ ニケーション能力はその伝達距離に応じて低下するものとする.ここで,骨梁表面の任意の点 xc までの情報伝 達 距 離 を l = xc − xb と し , 情 報 伝 達 可 能 な 領 域 Ω の 領 域 半 径 を lL と し た と き , l < lL で あ れ ば w(l ) = 3 (1 − l / lL ) π lL3 ,また l ≥ lL であれば w(l ) = 0 となる単調減少の重み関数 w(l ) を用いて,点 xc における力 学情報量 S sf を S sf ( x c ) = ∫Ω w(l ) Soc ( xb ) dΩ のように定義する.本来,骨細管の形態異方性を考慮すれば,細胞間ネットワークの情報伝達能力も異方性を 有すると考えられるが,本稿では,骨細胞は等方的な報伝達を仮定した.以上のように決定された力学刺激量 S sf ( xc ) をリモデリングの駆動力とし,この値が大きければ骨形成,小さければ骨吸収が行われるものとし,骨梁 表面移動速度 M& ( xc ) と力学刺激量 S sf ( xc ) との関係を Fig.7 に示すような正弦曲線で関連付けた. M& m ax = + d M& S sfC S sfL 0 Formation S sfU S sfB Resorption S sf M& m in = - d Fig.7: Relationship between mechanical signal Ssf and trabecular surface remodeling rate M& 5. 骨リモデリングシミュレーション 5.1 シミュレーションモデル 本報では,骨細胞の力学刺激感知機構と考えられている,細胞突起(骨細管)の異方性に着目し,等方的に 骨細管が分布しているモデル(Case A)と異方的に分布しているモデル(Case B)の比較を行った.シミュレーシ ョンモデルは,Fig.8 に示す単純な二次元骨梁モデルを採用した. Uniform displacement (2 MPa) 0.8mm Cellular process Cellular process Osteocyte Osteocyte trabecula trabecula surface X2 marrow X1 (Case A) trabecula surface (Case B) 0.8mm Fig.8: Single trabecular model, boundary condition, and cellular process morphology 解析領域(0.8 mm×0.8 mm 長方形空間)を一辺 8 μm の立方体 Voxel 要素で 100×100 に規則的に分割し た.また直径 160 μm の骨梁を荷重軸に対して斜めに配置し,空隙部は骨髄で満たされているものとした.骨 梁要素,および骨髄要素は等方線形弾性体と仮定し,材料定数は,それぞれ,Eb = 20 GPa,νb = 0.3,Em = 20 MPa,νm = 0.49 とした.また,Case A では,単位 Voxel あたり 10 本の骨細管が等方的に分布しているものとし, Case B では,10 本中2 本が骨梁表面に平行に,残りの 8 本は骨梁とは垂直な方向に分布していると仮定した. また,骨形成・骨吸収の役割を担う細胞が周囲の力学状態を感知する領域半径 lt は,240 μm とした.境界条 件として,下面の面外方向変位を拘束し,上面には 2 MPa の圧縮に相当する一様な変位を与えた.また,骨梁 の形態変化を表現する手法として,Level Set 法 9)を導入した. 5.2 シミュレーション結果および考察 シミュレーションにより得られた骨梁形態と圧力分布,圧力勾配方向の時間変化,および見かけの剛性,体 積分率の時間変化を Fig.9, 10 に示す.Case A, Case B いずれの場合においても,リモデリング初期段階にお いて,骨梁の鋭角部付近の圧力が局所的に高くなり,また曲げの影響を受けて,鈍角部付近の圧力が局所的 に低くなるため,この付近の圧力勾配が大きくなり,骨形成が優勢となった結果,骨梁全体的に太さを増した (10th step).その後,X1 方向中央部の骨梁が荷重支持の役割を担い,圧力が高くなる一方で,前述の鈍角部付 近の圧力は低下し,それに伴い圧力勾配も低下したため,骨吸収が優勢となった(20th step).その結果,骨梁 は最終的に荷重方向に配向した形態へと変化した. Initial state 10th step 20th step Case A 120th step Initial state 10th step 20th step 120th step Case B [MPa] Fig. 9: Change in morphology, pressure and its gradient vector Fig.10 に示すように,最終的な骨梁形態は,初期形態に比べ,骨梁の体積分率はほとんど変化していない か,あるいは減少しているのに対して,見かけの剛性は向上している.これは,骨梁形態の配向方向の変化に より,荷重支持機能が向上されたことを表しており,本研究で提案した骨リモデリングシミュレーション手法は, 実際に観察される骨の機能的適応現象を良く表現しているものと考えられる.また,Fig.10 に示すように,一連 のリモデリング活動を通して,骨梁内の大部分で圧力勾配ベクトルは骨梁の配向方向に対して垂直な成分が 大きい.その結果,Case A に比べ,Case B のほうが,骨細胞は力学刺激をより多く感知しており,見かけの剛 性や骨梁の体積分率においても,高い値を示している.これにより,骨細管の配向は,骨のリモデリング活動に 大きく影響を与えるパラメータの一つであると示唆される. 4.0 3.0 2.0 1.0 0.0 Case A Case B 0 20 40 60 step 80 100 120 Bone volume fraction, BV/TV Apparent stiffness σ2 / ε2 (GPa) 5.0 0.5 Case A Case B 0.45 0.4 0.35 0.3 0.25 0.2 0 20 40 60 step 80 100 120 Fig. 10: Apparent stiffness and bone volume fraction 6. おわりに 本研究では,骨小腔-骨細管系の構造と細胞レベルでの力学状態との関連性を明らかにすることを目的と した.まず,骨梁内の骨小腔・骨細管ネットワークの蛍光観察を行い,その形態的特徴について検討した.そ の結果,骨細胞間のネットワークは,骨梁表面に対して,垂直に形成されている様子が観察された.また,骨細 管の形態をファブリック楕円体を用いて評価した結果,ファブリック楕円体の主軸と骨梁表面に対する法線ベク トル報告がほぼ一致していたことから,骨細管の配向方向は,骨梁表面の法線方向に一致する傾向があること が示された.次に,骨細胞間のメカノセンサーネットワークの配向に着目し,ネットワークの形態的特長を考慮 した骨梁リモデリングの数理モデル化を行った.さらに,単純な二次元モデルを用いた骨梁リモデリングシミュ レーションを通じて,骨細管の配向が,骨リモデリング現象に与える影響について考察した.その結果,骨細管 を異方的に配置した骨梁モデルは,等方的なモデルに比べて,骨細胞が力学刺激をより多く感知しており,骨 細胞ネットワークの異方性が,骨のリモデリング活動に大きく影響を与える因子の一つであると示唆された. 本研究は,機械理工学専攻適応材料力学研究室の佐藤成道君の協力を得て進められた.ここに記して,謝 意を表す. 参考文献 (1) Sugawara, Y., Kamioka, H., Honjo, T., Tezuka, K., and Takano-Yamamoto, T., Three-dimensional reconstruction of chick calvarial osteocytes and their cell processes using confocal microscopy, Bone, 36, (2005), 877–883. (2) Knothe Tate, M. L., Adamson, J. R., Tami, A. E., and Bauer, T. W., The osteocyte, Int J Biochem & Cell Biol, 36-1, (2004), 1-8. (3) Parfitt, A. M., The Cellular basis of bone turnover and bone loss: A rebuttal of the osteocytic resorption: Bone flow theory, Clin Orthop & Related Res, 127, (1977), 236-247. (4) Cowin, S. C., Moss-Salentijn, L., and Moss, M. L., Candidates for the mechanosensory system in bone, J Biomech Eng, 113-2, (1991), 191-197. (5) 尾迫佑樹,安達泰治,田中基嗣,北條正樹,骨細胞によるメカノセンシングを考慮したリモデリングシミュ レーション,日本機械学会第 17 回バイオエンジニアリング講演会,04-48, (2005), 387-388. (6) 佐藤成道,安達泰治,田中基嗣,北條正樹,骨細胞ネットワークの形態的特徴を考慮した骨梁モデリング, 日本機械学会 2005 年度年次大会, 05-1-6, (2005), 13-14. (7) Adachi, T., Tsubota, K., Tomita, Y., and Hollister, S. J., Trabecular surface remodeling simulation for cancellous bone using microstructural voxel finite element models, J Biomech Eng, 123-5, (2001), 403-409. (8) Tsubota, K., Adachi, T., and Tomita, Y., Functional adaptation of cancellous bone in human proximal femur predicted by trabecular surface remodeling simulation toward uniform stress state, J Biomech, 35-12, (2002), 1541-1551. (9) Osher, S. and Sethian, J. A., “Fronts Propagating with Curvature- Dependent Speed: Algorithms Based on Hamilton-Jacobi Formulation”, J Comp Phys, 79, (1988), 12-49. !"#$%&'()*+,-./012345678 !"#$%&'(!")* +,-. Abstract: In the system-in-package (SiP) board, many devices made by several materials are embedded in a print circuit board, and the microstructures of print circuit boards cannot be ignored to evaluate the failure of embedded devices. In this study, the system to measure thermal deformation in micro region using digital image correlation in conjunction with an optical microscope was developed, and the system’s applicability to measurement of thermal deformation was verified by measuring purely thermal expansion of aluminum alloy. Based on the result, the distribution of strain in a SiP board during the thermal cycle test was measured using the developed system. Measured distribution of strain accurately related to the macro warpage of the circuit board. Furthermore, it was mentioned that the precision of measurement was affected with the image distortion, which often appears during these measurements, provided by the configuration of optical system Key words: Electronic Package, System in Package, Residual Stress, Digital Image Collation method 1.-9: -/0123456789:;<12=>?34@A34BC7DEF<GHIJ4KL8MN;OA PQRS&PQTUVWXOY8@AZ[\]^8_:;`\ab4cd<efgh_ijRkl8] :;mn[\opFRqrXs4tuvwxy`4zX{|y}_~<b4]^_BC@A3•€41 •4opF‚ijyc`ƒ`^„…X†v:;`\a‡ˆ‰w4opFŠ‹ijŒ•ƒy;4Ž••‘ ’“”•–(Digital Image Correlation Method: DICM)D<—˜™š›œX•ž8Ÿ _¡¢w_}<£" ¤4¥¦‚§¨\bƒw<©ª4«¬w4ijX-®ƒ_\¯°±²³´µ¶cd<·¸¹º_»¼½ 4q°X¾Fn[;`\a¿8Àfg˜Á‘<ÃĘÁ‘ƒ`:cB5ÅÆ4Š‹ij½4q°D ÇÈwÉ\[1]ayzy_Xn<b4]^8—˜™š4ÊË8•ž8›Ìu‚µ¶Í¼<ijÎu‚Ï ¦Ð¢\„ÑXÒÓ¦wÉ~<ÔÕZ[;`\²4w²Ö0.1pixels(1)~Ö0.01pixels[2]4]^8×Ø_ÙØ XÉ\ab4]^_ÒÚ¦_ DICM —˜™šXijŒ•ƒy;4Ûܱ‚m\cd8D<s4]^_ ÝÞßàšwáâXãäå\4zkæå\bƒXç„wÉ\a -è#$wD<àšéêàƒŽ••‘’“”•–‚ëFA|´c»–‚°`;<BC@A3‚Så\ 12=>?34opFŠ‹‚j¦å\ì–‚Ùãycaíy;<Ž••‘’“”•–—˜™š¿S4i jáâ4ãäÝÞßàš8¶`;îïy<í[n4ðñ4òó–ô]UÎuõJ8ö¢;4÷ø_s‚ ùúyca¿8£"¤8ûüycijáâ8¶`;D<ýþ˜™Â•‚°`;b4áâ‚ÿ!å\"x ì–‚Ùãy<í4S#±‚î$ycavc<è—˜™š4%°±‚î$å\cd8<1Ê&'4Àf gêÅÆw4opFijƒy;<()=>4*¼3RGHI?+=>4,3•€4opFŠ‹ij‚ô b_`<í[-[efg_·.ƒ4M"R FEM kl4/0ƒ4M"‚ôb_:ca 2.-!"#$;<%&'()*+456=7>?@A 2.1-!"#$;<%&'3BC- Ž••‘’“”•–D<Ž••‘’“J4©ª4’Ñ1823y< í41‚45ƒå\5Z_’“ÅÆ(678H9)?4:u;Š‹X<<4’“wDs4ÅÆ8Éc\z ‚:u;Š‹4”•8]~=>å\aíy;<b[n4?@‚M"å\bƒw2314§?‚m\ì– wÉ\aè#$wD<Sutton n[3]ƒA¹<678H9AB4”•4òó•C8D”D”•E‚=Fƒy c•C‚¦Gy<678H94§H8DIœ§HƒJ¹opFvwKLycMŽ‘‚N°ycab4§ HMŽ‘TUòó•CwD=>Ot8ô`;PQ_RS§CXT¶ü¸y<•CX•UHƒ_\cd< RS§C4Ϧ8D Newton-Raphson(N-R)–‚°`;`\ac¡y<N-R –8]~xy`k‚m\8D %V_WXY¦;XZ„ƒ_\a -íbw<678H94=>8]\§?ij4ì[ƒy;D<vp678H94\]ýþ4F‚KLy< Ÿ _òó•C‚°`c^=>‚ôb_^bƒw_’ÑŸ?4§?‚`aycaíy;b4/0‚²ƒ 8<Iœ§HTUJ¹opFvwKLycbC=>‚ôb_^bƒw<3h4_’Ñcd4†Îu_§ ?‚mca 2.2-DEF-678H94^=>OtwD<eâfg(Sequential Similarity Detection Algorithm: SSDA –)‚°`;_’ÑŸ?4§?‚m\ab4^=>Ot8ô`;<678H94§HD\]ýþ 4F‚KLy<678H94”•‚håòó•C8DE(1)4 S ;‚N°å\ab4 S ;DM"å\6 78H9iw4:u;4â‚håaS Xj5ƒ_\678H9‚WX678H9ƒAJ4²4ƒF_y< b[n4kl4â(u0, v0)‚23yc14§?ƒå\a -------- S ( u0 ,v 0 ) = M M " " I (X + u d 0 + i,Y + v 0 + j ) ! Iu ( X + i,Y + j ) -------(1) i= !M j=!M bbw<X, Y Dí[-[2314kl<i, j D231‚m1ƒå\678H9?nokl<M D678 H9p–<Iu<Id Dí[-[§Hqr4Ž••‘’“J4s’ÑXtå:u;‚håakl¥¦TU6 78H94`uv‚v 1 8tåa 2.3-G-EF-bC=w8ô`;<678H94”•4òó•C8D”D”•E‚²ƒ8¦GycE(2) 4 C ;‚N°å\a”•; C Dž8x4;‚x~<;X5Z`ys”•Xz`bƒ‚håa§Hr4 ’“Jw<C Xj5ƒ_\678H9‚WX678H9ƒAJ4678H9ƒF_y<§?{4§H| ‚m\a M ----- C ( u*,v *) = 1! M # # I ( X + u * +i,Y + v * + j) " I ( X + i,Y + j) d u ----(2) i= !M j=!M M M # # I ( X + u * +i,Y + v * + j ) d i= !M j=!M 2 " M M # # I ( X + i,Y + j) 2 u i= !M j=!M -678H9?4§?8DE(3)‚¦Gy<Iœ§HTUJ¹opFvwKLyc§HMŽ‘w678 H94=>‚ôb_^ab4ƒØ<(u, v)D678H94514§?‚hy<(u, v, !u/!x, !u/!y, !v/!x, !v/!y) D`p[²678H9?wJ¦z¶PQwÉ\cd<RSCD‡3wT¶ƒ_\a !u !u !v !v i+ j , v* = v + i + j -----------(3) ----------- u* = u + !x !y !x !y -òó•CX•UHwÉ\cd8<;Xj5ƒ_\678H94§H}~4=>8D N-R –‚°`< WX;8Dj5;•€¡ƒK¨n[\^=>4/0‚N°yca v 1 kl¤ƒ678H9MŽ‘4`uv 2.4-=7>?@AHI- Ž••‘’“”•–‚°`\bƒw<Jë4Ž••‘’“znòó1‚45 ƒå\678H94§H}~Xmn[\ayzy_Xn<678H94§HMŽ‘‚Ÿ 7y;`\c d<opF4•_s4ÎuX‚}<ij;‚ƒ@N°å\4D%„w_`aíbw<s678H945 14§?‚²ƒ8<ýþj5›…–8]~§?†‚•‡h·y<íbzn'œh¼8ô¢\©ª414 opF‚aˆycaýþj5›…–8]\•‡wD<S‰„Ñ–8]\?Šƒ‹_~<lè18Œ•Z [\bƒ_}<‡ÅÆw\Ž7Z[copF4•g•Š‹‚m\bƒXwØ\acJ<ij—˜™š4 ÊË‚v 2 8tåa 3JKLM3NO6PQRSTUVW 3.1 TU0X3BC--DICM 8]\§?ij8ô`;<1•4§?|ƒijZ[c§?|DZpy² v 2 ij—˜™š4`„ J‘y_`a’“4’XF4ðñw<'œh¼w4§?ƒ’“Jw4§?XJ‘y_`cdwÉ\aí bw<2314kl‚’“J4kl(xIm, yIm)zn'œh¼J4kl(xRe, yRe)½xy}§“y<ijZ[c §?|‚1•4§?|8§“å\bƒw"x‚ôb_^a -”âR•–"h—@4ðñ‚˜_HEwh·y<’XF4™x•ƒy;š¨c²4‚kl§“E(4) ƒy;¦Gyca "$ x = !x x id , x id = x Im + x dis ( x Im , y Im ) ----------- # Re ------------(4) $% y Re = !y y id , y id = y Im + y dis ( x Im , y Im ) bbw<(xdis, ydis)X’XF4™x•wÉ~<!x, !y Dij4«¬‚håa«¬•D©ª›œÝ•wÉ~< opF8ðñ‚•¨\4D(xid, yid)4FwÉ\a 3.2-VW3YZ -E(4)8ô`;<"x•Cƒ_\’XF4™x•(xdis, ydis)‚ž}8D<É\2314 ’“J4kl(xIm, yIm)ƒí[8Ÿqå\'œh¼J4kl(xRe, yRe)XÉ[ z`aíbw<¡¢£˜™Â •8]~¤S4§?|‚•¨<DICM 8]\ij;‚ƒM"å\bƒw™x•4¥¦ŽÂ•‚mca -¡¢£˜™Â•8]\\]ýþ§4§?ij‚ôb_`<mn[c§?‚’XF•C‚°`;"xy ca˜™Â•4ýþ|D("xRe, "yRe)=(-20µm, 20µm)ƒy<«¬D¨ 1.6(pixels/µm)wÉ:caijÅÆ 1400pixels©1400pixels ?4<"xqr4§?4Š‹v‚v 3 8tåa§?"x8]~’XF8]\ðñ Xÿ!Z[<ÅÆ?wJ¹_§?ƒ_:cavc<b4ƒØ4§?4 n¶ØD<(ªh_ijÎuw É\Šk®Ö0.02pixels 8”v\/0wÉ:c[3]a 4.-[01\3456=7 4.1-]^_`3a456=7- OAPQ‚«ö()=>‚Ÿ.ƒy;opFij‚ôb_:caOA PQ?3w4OY_Š‹X¬jZ[\cd<ij/04î$8Defg_¿±wÉ\-¬‚°`ca -()Ê&‚Så\<v 4 8tå/012°MŽ‘=>‚Ÿ.ƒy;<s®PQ4U¯°±C4âzn ä²\=>*¼4opFŠ‹‚ijycavp<ijw4WX}~‚Ϧå\cd<é³´µÂÞ˜ §?i‚°`;s¶u}~w4=>4Í~‚j¦yca-¬4¶u§7‚v 5 8tåab4ƒØ<-¬ XWX4¶u6·f‘ƒí[c¸46·f‘w‹_\¹º‚tå/0ƒ_:cab[D<¾»¼½Ë§ 4¾¿TUí[c¸4¹º8ÀÁ3Xðñ‚¢c/0wÉ~<¾»¼¶u8Ÿy;í4}~XJª8 DÏvn_`ayzy_Xn<JÃÀÁ4Äœ˜Åý1 TG cJ8šÆå\ƒ<y n}D¶u8Ÿy ;c¡J¶4-¬4¹º‚Ç\ac¸wD<ÈŸ4cd<šÆr4-¬¹º8É:;ijRKÊ‚ôb _:ca -v 5 Jw-¬XËgƒ_\¶u‚qr´Ì¶u TF ƒy<b4 TF ‚=Í}~ƒy;<h 1 8tå¶ u}~8ô`;ä²\opFŠ‹‚ DICM —˜™š8]~ijyca v 3-§?"xqr4Iœ§?8Ÿå\ij/0. v 4-™˜9=>4ÊË v 5-efg_=>4-¬4¶u§7 h 1-¶u¹º v 6-Captured image and measured region. -1•8ijwN°yc’“TUijÅÆ‚v 6 8tåaij«¬D 1.26(µm/pixels)wÉ:ca¶u¹ º B 8ô`;mn[copFŠ‹‚v 7 8tåa -΃opF#y D<ÏÐAB4@Ñ1wopFXÒ¨n[;`\²44<ÀÁ38ô`;ÏÐ8]\ y »õ4Œ•XÓ`cd<#x R# x 8MN;opF4;X×Ø`ab[n;45Z`#x R#xy DèÔ4OY_ Š‹8š¨;<áâ8]\ n¶ØX·|[;`\a -ÕÖP?34ÆopFD•ž8OY_Š‹‚ty;`\X<Äœ˜ÏÐ4•Xh_—@ƒA¹<op FŠ‹8²•X±X×n[\aíbw<ÏÐ4—@Xô]íJ•X8„\ÅÆ823y<?@ y 8ô¢ \ x »õ4΃opF#x ‚Ø„ÅÆw\Ùycab[8]~mn[copF!x 4 y Ú»õŠ‹Dv 8 8 É\]^8ƒU±‚tåab4bƒzn<efg8DÛ4Ÿ -ÜMŽ‘‚%°wØ\ƒK¨<v 8 ] ~mn[\•‡ƒU‚°`;=>4-¬‚aˆycas¶u¹ºwijZ[copFŠ‹‚²ƒ8y; mn[c-¬‚v 9 8tåab4;D<é³´µÂÞ˜§?i8]~efg8ijZ[;mn[c¬ƒJ‘y;`ca -JלêÝš8Š‹y;`\=>*¼4opFX<\Ùh8Defg_Í~ƒJ‘y;`\bƒzn< è—˜™šwijZ[cB5¾»¼J4opFŠ‹D<ÛÜ4ô¢\²4wÉ\ƒK¨n[\a 5.Jb: -è#$wD<Ž••‘’“”•–‚°`cBCÅÆ8ô¢\§?ÞopFij—˜™š‚Ùãyca -’“4’XF4]\áâDè#$wÙãyc"xì–‚°`\bƒw<ÿ!å\bƒXwØcab4 "xì–D½ßXÈà_¡¢w_}<"x•CX˜_HwÉ\cd8<s4]^_£"¤4ÊË8²% °wØ\ab4cd<"x‚ëFáÈ¡ DICM —˜™š8]\ijD<£"¤8ûüå\bƒ_}J ¦4Îu‚â$wØ\ƒK¨n[\a -cJ4ì–‚%V8ëFA|´\bƒ8]:;ijÎuÖ0.02pixels 4Šk®‚m\bƒXwØ\a -íy;<1Ê&'4ijwD<1•8ÀfgêÅÆw4ij‚ôb_`<OY_ÆopFŠ‹4ij X-®wÉ\zî$ycavp<()=>‚Ÿ.ƒycijwD<í4efg_Í~XJ‘y;`\b ƒzn<ij4ÛܱXâ$Z[\/0wÉ:ca v 7 State B 8Ÿå\ÆopF4ij/0 (a): ijãä, (b): #x., (c): #y, (d): #xy. ----v 8 Ùå7Z[copF4 y »õŠ‹- -v 9 Ùå7Z[copF]~…dn[\\Ù-¬ -----------------------ƒé³§?¤wijycefg_-¬4M" cdef [1] Ioannis Chasiotis, Mechanics of Thin Films and Microdevices, IEEE Transactions on Device and Materials Reliability, 4-2, (2004), pp.176-188. [2] Jinlong Chen, Ganmin Xia, Kebin Zhou, Guisuo Xia, Yuwen Qin, Optics and Lasers in Engineering, 43, (2005), pp.836-846. [3] H. A. Bruck, S. R. McNeill, M. A. Sutton and W. H. Peters III, Digital Image Correlation Using NewtonRaphson Method of Partial Differential Correlation, Experimental Mechanics, 29, (1989), pp.261-267. ------- 組成変調による薄膜材料の結晶構造制御 およびその電気・機械的特性の評価 工学研究科マイクロエンジニアリング専攻 神野 伊策 Abstract: We fabricate MEMS devices using functional thin film materials. The piezoelectric PZT films, which we have developed under control of their crystalline structure, are used to develop low voltage microactuators. In this study, piezoelectric MEMS switches are fabricated and the actuation properties are characterized. The introduction of the X-bar at the center of the beam causes the flat beam profile as well as the enhancement of the displacement at low voltage. Key words: MEMS, piezoelectric, PZT, thin films, RF-MEMS 1. はじめに 近年,MEMS 関連技術の研究の進展に伴い,各種マイクロアクチュエータやマイクロセンサ,およびこ れらを統合したデバイスの実用化が加速している.本研究では多様な機能性を有する強誘電体材料,そ の中でも特に圧電特性を用いた MEMS デバイス開発を目的として研究を行ってきた. 圧電材料は結晶構造 の異方性により応力による電荷の発生,ならびに電荷によるひずみの発生が生じる材料であり,電気機 械変換素子としてこれまで数多くの実用化がなされてきた.しかしながらその微細化・マイクロデバイ スとしての利用については,薄膜形成の困難さや信頼性の問題により実用例は限られたものになってい る.特に,多元素からなる複合酸化物を良好な結晶構造,一般にペロブスカイト構造を維持した状態で 薄膜成長させる必要があるが,基板拘束や格子のミスマッチ,また高温成膜に起因する熱応力等により その構造が大きくひずみ,特性が変化する.昨年度までに,PZT 圧電薄膜の特性制御を目的として,そ の結晶構造の詳細な評価を行った.はじめにエピタキシャル基板上に形成した c 軸配向エピタキシャル PZT 薄膜の結晶構造を 4 軸 X 線回折による基板上薄膜の結晶構造解析し,バルク材料の構造と比較した [1,2].更に基板から薄膜を剥離し,基板拘束をなくした状態での薄膜独自の構造を放射光 X 線による回 折および Rietveld 解析により評価し,薄膜の構造特異性は基板による影響のみではなく,プロセスに起 因した要素も存在することを明らかにした[3]. 本年は,薄膜材料の構造解析手法をベースに新しい圧電薄膜材料の開発,およびそれらを用いたデバ イス開発を行った.圧電薄膜材料として,今年度は相転移によるひずみの増大が期待できる反強誘電薄 膜の作製に取り組んだ. また非鉛材料として注目されている KNbO3-NaNbO3 系材料のエピタキシャル成長 をスパッタ法を用いて試み,配向膜の作製に成功した.このほか,PZT 圧電薄膜を用いた MEMS デバイス の開発を行い,特にミリ波帯域のスイッチングデバイスとして RF-MEMS スイッチ,および眼底検査の分 解能向上を目的とした収差補正用ディフォーマブルミラーの研究を行った.今回は特に応用デバイスと して汎用性の高い RF-MEMS スイッチ用圧電マイクロアクチュエータの研究について報告する. 高周波用スイッチングデバイスは,これまで半導体を用いたデバイスが主流であったが周波数が増加 しミリ波帯域にまで達するとデバイス内での損失が大きくなり,消費電力の増大等の問題が発生する. この対応としてメカニカルなスイッチング素子を MEMS 技術を用いて作製し, 挿入損失の改善を図る研究 が進められており, RF-MEMS スイッチとして注目されている[4]. 現在RF-MEMS スイッチの駆動源として, 静電力が主に用いられている.これは,電力消費量が小さい,高速レスポンス,微細加工が容易,とい った利点によるもので,実用化レベルの開発が進められている.しかしながら,静電引力によるスイッ チ駆動力は依然小さく, 駆動電圧が数十 V 必要となるため, 汎用の半導体素子による制御には向かない. 低電圧化の一つの取り組みとして,スイッチの駆動を圧電材料で行う方式が検討されている[5].圧電駆 動とすることで,大幅な駆動電圧の低下が可能となるが,これを実現するためには圧電特性の高い薄膜 (a) (b) Al Al Si Si 図1: Si 基板上に作製した PZT/Al 薄膜マイクロアクチュエータ: (a)片持ち梁構造, (b)両端支持梁構造 の形成が前提となる.今回,これまでの研究で 進めてきたRFスパッタ法を用いたPZT薄膜を用 いて圧電駆動 MEMS スイッチの試作を試みた. 2. スイッチの構造設計 スイッチの構造は圧電薄膜と振動板と呼ばれ る弾性層の2層からなるユニモルフアクチュエ ータ構造とした.圧電体の膜厚方向に電圧を印 図 2: 梁の中央部に X-bar 形状を有したマイクロアク 加することにより圧電横効果である電界に垂直 チュエータの構造 方向のひずみが発生するが振動板による拘束に より垂直方向の変位として振動する.はじめに アクチュエータ構造として片持ち梁および両端 支持梁構造のマイクロアクチュエータを試作し, その特性を調べた.図1にそれぞれの写真を示す. PZT薄膜の厚みは約2μmで,ほぼ同じ厚さのAl 層を電極兼振動板として形成した.長さ500μm の片持ち梁形状のマイクロアクチュエータの場 合,5Vで約1μmの変位が得られたが,各層の応 力バランスをとることが容易ではなく大きな初 期たわみが生じた.これを防ぐために長さ700 μmの両端支持梁構造のアクチュエータを作製 し,その写真を図1(b)に示す.図より初期たわ みのない良好な平坦性が確認できたが,変位量 は片持ち梁の場合の約1/10程度の低い値にとど 図 3: 1/2 モデルによるアクチュエータの FEM 変位 まり,剛性の増加による変位の低下が見られた. 解析結果 今回,片持ち梁構造の変位量を生かし,両端 支持梁の平坦性を両立させるアクチュエータ構 造について検討を行い,カンチレバー中央部にX形状の結合部を有する両端支持梁構造についてそのア クチュエータ特性を評価した.設計に際して,変位特性をMSC/MarcによるFEM解析を行い,その変位 量の評価を行った.今回設計したアクチュエータの構造を図2に示す.長さ350μm,幅200μmの2本の 片持ち梁を長さ100μmのX構造で結合した構造とすることで,膜中に堆積した残留応力の緩和およびそ り上がりを抑える効果を持たせた.アクチュエータ部の1/2モデルにおける計算結果を図3に示す.解析 の結果,5Vの電圧において約3μmの比較的大きな変位を発生することが可能であり,電送線路をアクチ ュエータの中心部を3μm以下のギャップで保持することにより5V程度の電圧でスイッチングが可能な デバイスが実現できることがわかる. 図4: Si 基板上に作製した X-bar 形状マイクロアクチ ュエータの外観写真 図5: 白色光干渉計により測定したマイクロアク チュエータの 3D 形状プロファイル 3. MEMS スイッチの試作 次に, 設計した構造を元にスイッチの試作を行った. 3.0 試作プロセスの概略は,①スパッタ法により Pt/Ti @ 0.1V 2.5 電極および PZT 薄膜を Si 基板上に成膜, ②振動板と (a) 2.0 して Cr 電極をスパッタで成膜, およびリフトオフに 1.5 よりパターンニング,③PZT/Pt/Ti をウエット,お よびドライエッチングにより梁構造にパターンニン 1.0 グ, ④裏面より Si 基板をドライエッチングしアクチ 0.5 ュエータの完成,以上のプロセスを経ることにより 0.0 3 4 4 4 圧電スイッチの試作を行った.完成したアクチュエ 1 10 1.5 10 2 10 0 5 10 Frequency (Hz) ータの写真を図 4 に示す.作製したアクチュエータ は単純な両端支持梁同様に平坦な形状を有しており, X-bar 構造による薄膜の内部応力解放および片持ち 2.0 梁のそり上がりを抑制することができた.詳細な形 (b) 状を白色光による干渉計により評価した.図 5 に測 1.5 定によって得られた 3 次元プロファイルを示す.測定結 果はz方向の形状を拡大して描画している.測定の結 1.0 果,完成した両端支持梁は支持部より最大で 3μm 支 0.5 持部より上に反った形状をしており,梁の幅方向にも最 大で 2μm の鞍型にたわんだ形状をしていることが明 0.0 0 5 10 15 20 らかとなった.初期たわみの絶対量は 3μm以下であ Volt (V) ることから,今後Cr振動板層の応力を制御することによ りその変形形状を更に減少させることが可能であると考 図6: X-bar 構造を有するマイクロアクチュエータの えられる. 変位特性: (a) 周波数依存性, (b) 電圧依存性 Displacement (µm) Displacement (µm) pp 4. アクチュエータ特性評価 引き続き作製したマイクロアクチュエータの変位特性を評価した.変位量はレーザドップラー振動計 を用いて行い,はじめに電極間に 0.1Vのサイン波を印加し,その周波数を 20kHz まで走引することで 共振周波数の測定を行った.図 6 に変位量の周波数依存性を示す.図より約 14.6kHz の周波数で共振が 確認できた.長さが 500μm の片持ち梁では共振周波数が約 7kHz であり,今回の設計でその共振は約 2 倍となった.片持ち梁の共振周波数は長さの 2 乗に反比例するが,今回の両端支持梁は長さが 350μm の片持ち梁が 2 つ組み合わさった構造をしていることから,共振周波数の増加は主に片持ち梁構造の長 さに起因したものと考えることができる. 次に共振周波数以下の 1kHz の周波数において,変位量と印加電圧との関係を調べた.結果を図 7 に 示す.図より変位量は電圧とともに増加し,15Vの電圧で約 1.5μmの変位が確認できた.しかしながら FEM の計算では 5V で約 3μmの変位量が得られていることから,今回の試作による結果は計算よりも 大幅に小さい値となった.理由として形成した PZT 薄膜の圧電特性が計算で用いた値(e31= -6.5C/m2) よりも小さい,また梁の形状が理想的な平面でないため変位に対する剛性が増加した,等の理由が考え ら得る.PZT 圧電薄膜の圧電特性を試作前にあらかじめ評価するとともに,振動板層の応力最適化によ り梁形状の平坦化を実現することで,理想的な振動が実現できると考えられる. 5. まとめ 本研究では,圧電薄膜を用いた MEMS デバイスとして RF-MEMS スイッチの開発を行い,設計,試作プロセ ス検討,およびアクチュエータ特性の評価を行った.Si 基板上に形成した PZT 薄膜を微細加工することで両端 支持梁構造のユニモルフマイクロアクチュエータを作製した.梁の中央部を X 形状とすることによって薄膜内 部の応力を緩和し,更に片持ち梁構造で見られた初期たわみを抑える効果が確認できた.形成したアクチュ エータの変位特性を評価した結果,共振周波数は 14.6kHz,変位量は 15V で約1.5μmの変位が確認できた. 今後,梁構造の最適化,および振動板層の応力を最適化することによって変位量の向上に取り組んでいく予 定である. 参考文献 [1] I. Kanno, H. Kotera, K. Wasa, T. Matsunaga, T. Kamada, R. Takayama, “Anomalous Crystalline Structure of Epitaxial Pb(Zr,Ti)O3 Films Grown on (100)Pt/(100)MgO by RF-Magnetron Sputtering”, J. Kore. Phys. Soc., 42 (2003) S1317 [2] I. Kanno, H. Kotera, K. Wasa, T. Matsunaga, T. Kamada., R. Takayama, “Crystallographic characterization of epitaxial Pb(Zr,Ti)O3 films with different Zr/Ti ratio grown by radio-frequency-magnetron sputtering”, J. Appl. Phys. 93, (2003)4091 [3] I. Kanno, H. Kotera, T. Matsunaga, K. Wasa, “Intrinsic crystalline structure of epitaxial Pb(Zr,Ti)O3 thin films”, J. Appl. Phys., J. Appl. Phys., 97, (2005) 074101 [4] G. M. Rebeiz, “RF MEMS theory, design, and technology”, (2003), Wiley-Interscience [5] H. Lee, J. Park, J. Bu, “Piezoelectrically Actuated RF MEMS DC Contact Switches With Low Voltage Operation”, IEEE Microw. Wireless Compon. Lett., 15 (2005) pp.202-204 誘電体ナノスプリングを有するコンデンサ型デバイスの開発 とナノ電気機械要素への応用 工学研究科マイクロエンジニアリング専攻 鈴木 基史 Abstract: We have demonstrated the capacitors with an insulator layer which consists of nanosprings created by dynamic oblique deposition. The capacitors show the rather large leakage current when they are in air. However, the onset voltage of the current-increase increases with repeating the voltage scans. Furthermore, the leakage current decreases significantly under the vacuum higher than ½¼ Torr. This indicates that the leakage current is due to the gas and/or humidity absorbed on the porous nanospring layers. We tried to apply these capacitors to thin film actuators and humidity sensors. The leakage current under the constant bias voltage increases with increasing humidity, although the separation between the electrodes did not change. Therefore, the capacitors with porous insulator are possible candidate of the noble humidity sensors, while the much effort may be required to realize the thin film actuators Key words: dynamic oblique deposition, thin films, nanotechnology 1. はじめに 我々は,真空蒸着による成膜のプロセスにおいて,基板と蒸着源の方位関係を制御する動的斜め 蒸着によって薄膜の内部形態を制御し,形態と密接に関わる機能性の発現と,その応用に関する 研究を進めている.動的斜め蒸着法によれば,成膜中の基板の面内方位の制御によって螺旋やジ グザグなどの様々な形態を成長させることが可能でるが [1],共通したかたちの特長は「細長い」 ことであり,これに起因する光学的,磁気的な異方性を示すことが知られている.当然ながら機 械的な性質に対しても興味がもたれているが,これまでの研究例ではバネ形のコラム構造の膜厚 「細長い」かたちで最も興味深い横 方向への準静的な伸縮特性が評価されているのみである [2, 3]. 方向の変形や,能動的な応用に関する研究はほとんど進んでいない. そこで我々は本助成研究で動的斜め蒸着膜内部に形成されるナノコラムの横方向の変形特性の 測定や,外部信号に応答して変形する新しいアクチュエータの実現を目的とした研究を開始した. 昨年度までにパターニング基板上へのナノコラムの選択成長に成功した.本年度はこれを,本拠 点において微小構造体の横方向の変形特性を評価するための独自技術を有する物性工学講座・材 料物性学分野 (北村研究室) に引き継ぎ,世界に先駆けてナノコラムの横方向の変形特性の評価へ と展開した. 一方我々は,研究開始当初からの「夢」である外部からの信号によって「動く薄膜」の実現に向 けて,螺旋型のコラム構造 (ナノスプリング) を有する高耐圧コンデンサの開発に着手した.その 結果,静電気力による薄膜の伸縮実現にあと一歩にせまる高耐圧を実現した.また,検討の途中 でコンデンサの電極間を流れる電流が,周囲の環境の湿度に敏感に応答するという思わぬ副産物 を得た.本報告書では,圧電材料を一切用いない新奇な薄膜アクチュエータのアイディアと,本 年度実施した実験結果を報告する. 2. 薄膜アクチュエータのアイディアと目標の設定 図 1 に示すように,2 枚の金属で誘電体をサンドイッチした平行平板コンデンサにおいて,誘電 体層がナノスプリングの集合体で構成されていれば,電極間に電圧を加えることによって発生す る静電引力によって,膜厚を縮めることができるはずである.我々はこのような積層構造を動的 斜め蒸着法によって作製し,圧電材料を用いない新しいタイプのアクチュエータの実現を目指し ている.本節では,そのために必要な膜の特性について考察する. ナノスプリングのばね定数を とすると,ばねを だけ変形させるために必要な力は で ある.一方,一本のナノスプリングが占める面積を , 長さを ,有効比誘電率を とすると,電 極 Al電 キャップ層 電極 ナノスプリング層 誘電体ナノスプリング 有効比誘電率: ε 厚さ: d (Ta2O5) シード層 ITO電極 ガラ ス基 板 電極 図 1: 誘電体ナノスプリングを有するコ ンデンサの模式図. 圧 を加えたときに発生する力は 図 2: 誘電体ナノスプリングを有するコンデンサの SEM 像. であるから, (1) ¢ である.過去の動的斜め蒸着膜の研究から, m, nm , N/m,また誘電 とすれば,50 V 印加時 体に比較的大きな比誘電率を有する Ta O を用いることを想定して の変位量は 2.6 nm である.この値は原子間力顕微鏡や,偏光解析法などの光学的手法によって十 分検出可能である.一般の成膜法によって成膜された高密度で均一な酸化物薄膜であれば,この 程度の耐圧を実現することは不可能ではないが,動的斜め蒸着膜は低密度で不均一な内部構造を 有するため,十分な絶縁耐圧を確保できない可能性がある.そもそも動的斜め蒸着膜を誘電体層 に用いたコンデンサ自体,これまでほとんど作製された例がなく,電圧印加時の電流-電圧特性な ど,基本的な特性が全くわかってない.そこで本研究では,電極間距離が 1 m 程度で,50 V の 電圧で絶縁破壊しない誘電体ナノスプリングを有するコンデンサ型デバイスを開発することを第 1 の目標にした. 3. ナノスプリングを有するコンデンサの作製 本研究では最終的にコンデンサの変位を AFM または偏光解析法で測定することを想定し,電極 の一方に透明導電膜を採用することにした.ガラス板上に 3 mm 幅の ITO 薄膜を形成した基板を用 意し,その上に Ta O ナノスプリングを形成した.平坦な基板上に多孔質のナノスプリング層を 直接形成すると,電極と誘電体層の界面の凸部に電界が集中して絶縁破壊しやすくなることを懸 念し,ナノスプリングの上下を高密度の Ta O 層でサンドイッチすることにした.図 2 に典型的 なコンデンサの断面の SEM 写真を示す.まず基板を蒸発源に対して正対させて 300 nm の Ta O 蒸着し,緻密なシード層を形成した.続いて蒸着角を Æ に保持し,蒸着速度に比例した回転速度 で基板をゆっくり面内回転させながら成膜した.回転速度は形状膜厚が 850 nm 成膜される毎に 1 回転の割合とした.図 2 に示した試料ではナノスプリングを形成する間に,面内角を 1.5 回転 回転 させて作製した.太さ 100 nn 程度の細い長いコラムが途切れることなく螺旋を形成していること がわかる.このナノスプリング層の上に,再び基板垂直方向から Ta O を蒸着し,高密度のキャッ プ層を形成した.ナノスプリングから均一層への遷移に必要な厚さは 130 nm で,キャップ層はこ れよりも十分厚く 300 nm の厚さにした.最後にメカニカルマスクを用いて 3 mm 幅の Al 電極を 形成した.ナノスプリング層は機械的に非常に弱いため,メカニカルマスクは試料に絶対に接触 しないようマスクホルダを工夫した. 5.0x10-5 4.0x10-5 6.0x10-7 air vacuum -7 5.0x10 -7 4.0x10 current(A) current(A) -5 3.0x10 -5 2.0x10 1.0x10-5 2.0x10-7 0.0x100 -1.0x10-5 -80 3.0x10-7 1.0x10-7 -60 -40 -20 0 20 40 60 80 0 0.0x10 20 30 40 50 voltage(V) 図 3: 誘電体ナノスプリングを有するコンデン サの I–V 特性.真空中と大気中での測定結果. 60 70 80 90 RH(%) 図 4: 相対湿度と素子を流れる電流の関係 4. ナノスプリングを有するコンデンサの電流–電圧特性 ナノスプリングを有するコンデンサの電流–電圧特性や,絶縁耐圧などは研究報告がなく全く未 知である.我々がアクチュエータを実現するためには,まずこれらを明らかにする必要がある.図 3 は,図 2 に示した試料の ITO 側の電位を横軸に,縦軸にコンデンサに流れる電流値をとった電 流–電圧特性である.大気中に素子をおいて電圧を徐々に増加させると,はじめて 10 V 前後に達し たときに突然大きな電流が流れはじめる.しかしながら電圧を下げれば若干の履歴は示すものの ほぼ可逆的に電流が減少するため,絶縁破壊ではない.同様の電流増加が負バイアス側でもおき る.電圧を正負に掃引しながら振幅を次第に大きくしていくと,電流の立上り電圧が次第に大き くなることがわかった.多孔質のナノスプリングに吸着したガスや水分が電流を担っている可能 Torr 台の真空中で電流– 性があるため,これを検証するためにコンデンサを真空槽にいれて 電圧特性を測定すると,図 3 の実線で示したようにバイアスが 60 V を越えても大きな電流が流れ ないことがわかった. 以上 I–V 測定の結果から,少なくとも真空中であれば十分な耐圧を達成できるし,空気中でも 電圧を何度か正負に掃引すれば目標の耐圧に近い高耐圧が実現できることがわかった.また,空 気中で流れる大きな電流は,空気や水分などによるものであることはまちがいなく,この特性を 応用したガスセンサや湿度センサ実現できる可能性を示唆している.次節において,これらの応 用検討した結果を簡潔に述べる. 5. ナノスプリングを有するコンデンサの応用 5.1 アクチュエータ 目標とした耐圧を実現できたので,AFM と偏光解析法を用いて電圧印加 時の電極間距離の変化を調べた.測定は装置の制約から全て大気中で行った.両手法とも先に想 定した数 nm の変化があれば十分検出可能な方法であるが,残念ながら本年度作製した試料では, 電圧印加による有意な変形は検出することができなかった.その原因としては,作製したナノス プリングの密度が高く,隣り合うスプリングの干渉が変形の妨げとなったこと,ナノスプリング は反電場の効果によって長手方向に分極しやすく,そこに働くトルクが膜厚の収縮を打ち消す方 向に働いたことなどが考えられる. 5.2 湿度センサ ナノスプリングに吸着したガスや水分が電流に寄与することを利用し,湿度セ ンサへの応用の可能性を検討した.図 4 はコンデンサに一定バイアスを印加して相対湿度を変化 させた際の電流の変化を示す.湿度の増加とともに電流が単調に増加することがわかる.また図 には相対湿度を増加させるとき () と減少させるとき (¯) で大きな差がなく,湿度センサとして有 望であることがわかった. 6. まとめと今後の展開 誘電体ナノスプリングを 2 枚の電極でサンドイッチしたコンデンサを作製することに成功した. 電流–電圧特性を測定したところ,大気中ではガスや水分の吸着に起因すると思われる電流が顕著 であったが,バイアスを正負に何度か掃引するうちに,50 V 程度の大きな耐圧が得られることが わかった.一方真空中では 80 V 以上の大きなバイアスを印加してもほとんど電流が流れず,大気 中に比べて大きな耐圧が実現できることがわかった. 大気中で電圧を印加した際の膜厚の変化を調べたが,残念ながら有為な膜厚の変化を検出する ことはできなかった.しかしながらこのコンデンサを流れる電流が湿度に敏感に応答することを 見出し,新奇湿度センサという思わぬ副産物を得た.今後湿度センサについては時期をみて実用 化の検討を行う予定である.また,アクチュエータについては本年度の検討でうまくいかなかっ た原因をさらによく吟味し,装置の改良などを含めて進め方を再検討したいと考える. 謝辞 本研究は岸本 一昭,中嶋 薫 助手,木村健二 教授の協力を得て推進した.また,本学マイクロ エンジニアリング専攻の小寺研究室の皆様,島研究室の皆様,木下技官,機械理工学専攻の北村 研究室の皆様には大変お世話になりました.ここに御礼申し上げます. 研究発表実績 (関連の研究成果を含む) 口頭発表: 国内会議: 5 件,国際会議: 5 件 学術論文, 会議論文: ¯ M. Suzuki, W. Maekita, Y. Wada, K. Nakajima, K. Kimura, T. Fukuoka and Y. Mori, “In-line aligned and bottom-up Ag nanorods for surface-enhanced Raman spectroscopy,” Appl. Phys. Lett., in press. ¯ M. Suzuki, Y. Wada, W. Maekita, K. Nakajima, K. Kimura, T. Fukuoka and Y. Mori, “Physically SelfAssembled Ag Nanorod Arrays for Tunable Plasmonic Sensors,” e-J. Surf. Sci. Nanotech., 3 (2005) 280283. ¯ M. Suzuki, W. Maekita, Y. Wada, K. Nakajima, K. Kimura, T. Fukuoka and Y. Mori, “Surface Enhanced Raman Scattering on Physically Self-Assembled Ag Nanorod Arrays,” Mater. Res. Soc. Symp. Proc. 900E (2006) 0900-O08-03. ¯ M. Suzuki, W. Maekita, K. Nakajima, K. Kimura, T. Fukuoka, Y. Mori, “Tunable Plasma Resonance of Ag Nanorod Arrays Prepared by Dynamic Oblique Deposition,” 2005 MRS Fall Meeting, AbstractViewer (2005) Abstract No. GG3.2. ¯ T. Fukuoka, M. Suzuki, Y. Wada, K. Nakajima, K. Kimura, Y. Mori, “Surface-enhanced Raman Scattering of Rhodamine 6G Using Physically Self-assembled Ag Nanorod Arrays,” 2005 MRS Fall Meeting, AbstractViewer (2005) Abstract No. Ra5.10/Rb5.10. ¯ T. Fukuoka, R. Kuramoto, N. Shinohara, S. Ikuta, M. Suzuki, Y. Mori, “Stable and Reproducible Chainlike Assembly of Gold Nanoparticles for Surface-enhanced Raman Scattering,” 2005 MRS Fall Meeting, AbstractViewer (2005) Abstract No. Abstract No. O12.30. 特許出願: 1 件 (JST の支援を受けての PCT 出願) 参考文献 [1] M. Suzuki, Y. Taga, Integrated sculptured thin films, Jpn. J. Appl. Phys. Part 2 40-4A (2001) L358– L359. [2] M. Seto, K. Robbie, D. Vick, M. Brett, L. Kuhn, Mechanical response of thin films with helical microstructures, J. Vac. Sci. & Technol. B 17-5 (1999) 2172–2177. [3] D. Liu, D. Ye, F. Khan, F. Tang, B. Lim, R. Picu, G. Wang, T. Lu, Mechanics of patterned helical Si springs on Si substrate, J. Nanoscience and Nanotechnology 3-6 (2003) 492–495. 生体組織マトリクスの変形・損傷場と細胞ネットワークの相互作用 工学研究科機械理工学専攻 田中 基嗣 Abstract: The osteocyte network system embedded in bone matrix is considered as one of the important elements in the mechanosensory mechanism of the bone remodeling. This study observed the calcium response of isolated osteocytes and osteocytes embedded in the living bone tissue matrix to the direct deformation by the glass micro needle. To compare the short-time mechanosensory response of osteocytes and osteoblasts, we observed the calcium response also in isolated osteoblasts to the direct deformation. The changes in the intracellular Ca2+ concentration were observed using a confocal laser-scanning microscope. As a result, it was shown that the osteocyte cell body was less sensitive than the cell process. In addition, the role of cytoskeletal actin filaments in the mechanosensory mechanism was discussed in terms of its localized distribution. It was suggested that the mechanosensory mechanism in osteocytes is related with the cytoskeletal actin filament structure. Finally, we could observe the calcium response of osteocytes embedded in the living bone tissue matrix in situ. This knowledge can be applied to clarify the relation between the macroscopic mechanical stimulation to the bone matrix and the microscopic mechanosensing of osteocytes in detail. Key words: Osteocyte network, Mechanosensory mechanism, Living bone tissue matrix, Direct deformation 1. はじめに 骨細胞は,骨組織表面において骨形成を担う骨芽細胞が自ら分泌した骨基質中に埋め込まれて分化し たものであり[1],ギャップジャンクションによりお互いに結合することで骨基質内で複雑な骨細胞ネッ トワークを形成している[2].骨細胞は,骨梁に対する力学負荷により生じる骨基質の変形・損傷場を感 知し,骨細胞ネットワークを介して刺激情報を周囲の細胞に伝達する機能を有すると考えられている[3]. 実際に,in vivo 実験によって,骨細胞が力学負荷を感知する細胞であることが確認されている[4].硬い 骨基質中に埋め込まれて存在する骨細胞に関する研究は,他の骨系細胞と比べて十分ではなかったが, 骨細胞単離法の確立[5]および骨細胞樹立株の確立[6]により,生化学・力学刺激に対する細胞応答の観察 例[7-11]が報告されるようになった.特に,細胞全体に対して一様に与えられた生化学・力学刺激に対し て,カルシウム応答が発現することが報告されている[8, 9].しかしながら,従来は,細胞全体に一様に 与えられた刺激に対する細胞応答の検討が主であり,単一の骨細胞における刺激付与部位と細胞応答を 詳細に関連付けた報告例は少ない.そのため,骨細胞の詳細な力学刺激感知・情報伝達機構は解明され ていないのが現状である.また,骨細胞樹立株 MLO-Y4 細胞が細胞体内に明瞭なアクチン骨格構造を有 する[12]のに対し,生体から単離した骨細胞においては細胞突起のみにアクチン骨格構造が局在する[13] ことが知られている.多くの細胞において,アクチン骨格構造は,力学刺激感知機構に重要な役割を果 たすことが示唆されている[14]ことを考えると,単離骨細胞を用いた研究が望まれる.また,生体内で 骨細胞は周囲を硬い骨基質で覆われているため,骨のリモデリング現象における力学的階層構造を考え ると,生骨組織へのマクロな力学刺激とミクロな骨細胞の刺激感知を結びつけた研究が望まれる. そこで,本研究では,骨細胞の複雑な力学刺激感知メカニズムの解明に向けた第一歩として,骨組織 から単離した骨細胞における刺激付与部位と細胞応答を関連付けること,および生骨組織へのマクロな 力学刺激に対するミクロな細胞応答を観察することを目的とした.まず,骨組織から単離・培養した骨 細胞に対してマイクロニードルにより局所変形を与え,その際の細胞応答としての細胞内 Ca2+濃度変化 をその場観察し,単離骨細胞の力学刺激応答と細胞内アクチン骨格構造との関連について検討した [15-18].その結果得られた骨細胞における応答特性を,分化前の形態である骨芽細胞(骨組織から単離・ 培養したもの)の応答特性と比較した.また,生骨組織においてマイクロニードルによる骨基質への局 所変形を与えた際に生じる骨細胞のカルシウム応答を観察した[19-20]. 2. 実験方法 2.1 骨組織採取・培養方法 本研究では,13日齢ニワトリ胚より採取した頭蓋冠を用いた.骨組織表面 の骨膜中には,骨芽細胞をはじめとする骨細胞以外の細胞が存在するため,1 mg/mlのCollagenase Type I (SIGMA) 溶液中で30分振盪培養(37.5°C環境下)し,骨膜を除去した.次に,0.05 % EDTA (Gibco) 溶 液中で15分間振盪培養(37.5°C環境下)する脱灰処理を2回施した.以上のようにして得られた生骨組織 を,α-MEM (Gibco) に2 % FBS (SIGMA) および抗生物質を添加した培養液中において,12時間培養(温 度37°C,湿度100 %,5 % CO2-95%Air環境下)した. 2.2 骨細胞単離・培養方法 骨細胞単離・培養方法 13 日齢ニワトリ胚より採取した頭蓋冠を用いて, Kamioka ら[13]の手法に準じ,骨細胞を単離した.単離した骨細胞を,poly-D-lysine(SIGMA)および fibronectin(COSMO BIO)でコーティングしたガラスボトムディッシュ(φ = 35 mm,MatTek)に,5.0x103 ~1.0x104 cells/dish の細胞密度で播種した.培地として 10% FBS 溶液を加えたα-MEM を用い,15~18 時間培養(温度 37°C,湿度 100%,5%CO2-95%Air 環境下)した. 2.3 骨芽細胞単離・培養方法 本研究では,骨細胞のアクチン構造骨格と細胞応答の相関を調べるにあ たって,分化前の形態である骨芽細胞についても同様の検討を行い,両者の応答特性を比較した.2.1 節で採取した骨組織表層の骨膜中に存在する骨芽細胞を単離し,2.2 節と同様にして培養した. 2.4 単離骨細胞の力学刺激応答観察方法 単一の単離骨細胞に対して,先端部を直径約 2 μm の滑らか な状態に加工したガラスマイクロニードルを直接押し込んでスライドさせることにより,局所変形を与 えた(図 1(a)) .多くの細胞において力学刺激感知と深い関連性が示唆されている[14]アクチンフィラメ ントに着目し,その分布と細胞応答の相関を調べた.単離骨細胞においては,アクチンフィラメントが ほとんど存在しない細胞体部分とアクチンフィラメントが局在する細胞突起(根元)部分に対して局所 変形を与えた.局所変形に対する細胞応答としての細胞内 Ca2+濃度変化を観察するため,[Ca2+]i 蛍光指 示薬 Fluo-4 AM(1 mM solution in DMSO, Molecular Probes)を細胞内に導入した.観察には,共焦点レーザ ー走査型蛍光顕微鏡(LSM510, Zeiss)を用い,100 倍の油浸対物レンズを用いた.単離骨芽細胞につい ても,同様にして,マイクロニードルによる局所変形に対するカルシウム応答を観察した.骨芽細胞に は細胞突起が存在せず,アクチンフィラメントが細胞体内部に分布・配向して存在することが知られて いる[21]ため,単離骨芽細胞に対しては,細胞体に局所変形を与えた. 2.5 生骨組織中の骨細胞の力学刺激応答観察方法 カバーガラス上に固定した生骨組織に対して,先端 直径約 2 μm のガラスマイクロニードルを押し込 むことにより,骨細胞周辺の組織に局所変形を与 えた(図 1(b)) .観察面から約 45°の角度でマニピ Nucleus Osteocyte, ュレータ (NI2, eppendorf) に取り付けたマイクロ Osteoblast 1~2μm ニードルをその軸方向に動かすことにより,組織 中の骨細胞に近づけた.この際,一種類の蛍光指 (a) 単離骨細胞・骨芽細胞への局所変形付与 示薬のみを用いると,局所変形付与時に細胞の高 さの変化やフォーカスのずれによる蛍光輝度の変 化が生じ細胞応答観察が困難となるため, Ratiometry 法[22]を用いた.Ratiometry 法は,Ca2+ 濃度を二種類の異なる蛍光指示薬の蛍光輝度比で 表す手法である.これにより,蛍光輝度観察に与 える影響を低減できる.細胞内 Ca2+濃度が上昇す ると蛍光輝度が上昇する Fluo 4 と低下する Fura Red を細胞内に導入し,これらの蛍光輝度比を測 定した.観察には,共焦点レーザー走査型蛍光顕 微鏡 (LSM510, Carl Zeiss) を用い,100 倍の油浸対 (b) 生骨組織への局所変形付与 物レンズを用いた. 図 1 局所変形付与方法 (c) (b) (a) 10μm 10μm 10μm (a) t = -5 sec (b) t = 0 sec (c) t = 10 sec 図 2 単離骨細胞細胞体への局所変形付与にともなう細胞内 Ca2+濃度変化 (a) (b) 10μm (c) 10μm 10μm 10 40 5 30 0 20 10 Mechanical stimulation -5 0 5 10 15 Time, t (sec) 20 60 15 50 -5 -10 25 図 4 単離骨細胞細胞体への局所変形付与におけ る蛍光輝度およびマイクロニードル先端位 置の時間変化 Fluorescent intensity Sliding distance of microneedle Peak of fluorescent intensity Sliding distance of microneedle, δ (mm) Fluorescent intensity, F 60 Fluorescent intensity, F Fluorescent intensity Sliding distance of microneedle Sliding distance of microneedle, δ (mm) (a) t = -3 sec (b) t = 0 sec (c) t = 10 sec 図 3 単離骨細胞細胞突起への局所変形付与にともなう細胞内 Ca2+濃度変化 20 50 15 40 10 30 5 20 10 0 Mechanical stimulation 0 5 10 15 Time, t (sec) 20 -5 図 5 単離骨細胞細胞突起への局所変形付与にお ける蛍光輝度およびマイクロニードル先 端位置の時間変化 3. 結果および考察 3.1 単離骨細胞および単離骨芽細胞の力学刺激応答 単離骨細胞に局所変形を与えた前後の蛍光輝度 観察写真を,図 2(細胞体への刺激)および図 3(細胞突起根元部分への刺激)に示す.また,図 2(a) および図 3(a)中に四角で囲んだ領域内の平均蛍光輝度の時間変化を,図 4(細胞体への刺激)および図 5 (細胞突起への刺激)に示す.実線は蛍光輝度を示し,破線は透過光観察写真から測定したマイクロニ ードル先端位置を示す.なお,マイクロニードルで細胞に触れた瞬間を時刻 t = 0 sec およびマイクロニ ードル位置 δ = 0 mm とした.細胞体部分に局所変形を与えた場合,変形付与直後の急激な蛍光輝度変 化は観察されなかった(図 2(b)) .その後,マイクロニードル押込み量に対応して,蛍光輝度のゆるやか な上昇・下降が見られた(図 2(c)) .この蛍光輝度の上昇・下降は,細胞応答としての細胞内 Ca2+濃度変 化を示すものではなく,マイクロニードル押込みによる細胞体高さの変化によって生じる高さ方向の映 (a) 20μm (b) 20μm (c) 20μm (a) t = -3 sec (b) t = 0 sec (c) t = 10 sec 図 6 単離骨芽細胞への局所変形付与にともなう細胞内 Ca2+濃度変化 Sliding distance of microneedle, δ (mm) Fluorescent intensity, F Fluorescent intensity り込みが主原因であると考えられる.以上より,単離骨 Sliding distance of microneedle 細胞においては,細胞体への局所変形に対して,顕著な 80 20 Peak of fluorescent カルシウム応答が発現しないことがわかった. 15 70 intensity 次に,単離骨細胞の細胞突起根元部分に局所変形を付 10 60 与した場合,マイクロニードルが細胞突起根元部分に触 5 れたと同時に蛍光輝度が急激に上昇した(図 3(b)) .細胞 50 0 体高さの変化なしに,急激な輝度上昇が見られたことか -5 40 ら,単離骨細胞においては,細胞突起根元部分の局所変 -10 Mechanical 30 stimulation 形に対してカルシウム応答が発現することがわかった. -15 その後,引き続いてマイクロニードルを細胞体部分にま 20 -20 -10 -5 0 5 10 15 で押込むと,マイクロニードルの押込み量に対応した蛍 Time, t (sec) 光輝度の上昇・下降が見られた(図 3(c)) .この蛍光輝度 図 7 単離骨芽細胞への局所変形付与にお の上昇・下降は,マイクロニードル押込みによる細胞体 ける蛍光輝度およびマイクロニード 高さの変化によって生じる高さ方向の映り込みが主原因 ル先端位置の時間変化 であると考えられる.マイクロニードルを単離骨細胞か ら離した後も,蛍光輝度は,変形付与前のレベルまで下降しなかった(図 5) . 比較の対象として,単離骨芽細胞の局所変形付与前後の蛍光輝度観察写真を図 6 に,図中四角で囲ん だ領域内の平均蛍光輝度の時間変化を図 7 に示す.細胞体への局所変形に対する単離骨細胞の応答特性 とは異なり,単離骨芽細胞は,細胞体への局所変形に対して急激なカルシウム応答を示した. 3.2 アクチン骨格構造と細胞応答の相関 3.1 節より,細胞体内部にアクチンフィラメントが配向して 存在する単離骨芽細胞は,細胞体に対する局所変形に対して敏感にカルシウム応答を発現させることが わかった.一方,アクチンフィラメントが細胞突起に局在する単離骨細胞においては,細胞体への局所 変形に対する顕著なカルシウム応答は発現せず,細胞突起根元部分への局所変形に対してカルシウム応 答が発現することがわかった.以上から,単離骨細胞は,アクチンフィラメントが局在する細胞突起へ の局所変形に対する応答特性が高いことが示唆される. また,アクチンフィラメントは,細胞に付与された力学刺激の伝達経路としての働きを担っていると 考えられている[14].骨芽細胞は骨梁表面に存在するために,細胞体内部にアクチンフィラメントを配 向させ,骨梁表面に強く接着する必要があると考えられる.一方,骨細胞は,間質液で満たされた骨小 腔内に存在する. そのため, 細胞体内部にアクチンフィラメントを配向させる必要がないと考えられる. また,骨細胞は,隣接する骨細胞と,骨細管内に伸張した細胞突起同士をギャップ結合を介して結合し 合っている.そのため,細胞突起内部には,アクチンフィラメントを局在させる必要があると考えられ る.これらのことより,骨芽細胞が自ら形成した骨基質中に埋没する過程で骨細胞へと分化する際に, 周囲の力学環境の変化にともなって,その環境下で細胞形状および細胞間ネットワークを構築・維持す るためにアクチン骨格構造を適応的に変化させることが示唆される.これにより,骨芽細胞と骨細胞の 応答特性に違いが生じることが考えられる. Fluorescence ratio, R 3.3 生骨組織中の骨細胞の力学刺激応答観察 生骨組織に対して骨細胞周辺部位に局所変形 を加えた際の骨細胞の蛍光輝度比変化観察結果 を図 8 に示す.図中では,蛍光輝度比を擬似的 に 256 段階の輝度に変換して表示した.また, 5μm 図中四角で囲んだ領域内における平均蛍光輝度 比の時間変化を図 9 に示す.ニードル先端が組 (a) t=-10 sec (b) t=0 sec (c) t=10 sec 織内の標的骨細胞に最接近した時点を,時刻 t = 図 8 生骨組織への局所変形付与にともなう骨細胞 0 sec とした.本結果より,生骨組織への局所変 内 Ca2+濃度変化 2 形に対して,生組織内の骨細胞のカルシウム応 答を初めて確認できた.これにより,骨細胞近 1.5 傍の骨基質の変形と細胞応答の相関を検討する 上で, 新たな知見を得ることができたと言える. 1 また,単離骨細胞においてアクチンフィラメ ントとカルシウム応答の深い関連性が示唆され Closest approach 0.5 たこと(3.2 節)を考慮すると,今後,生骨組織 of microneedle 中の骨細胞のカルシウム応答と同時に,細胞内 0 -20 0 20 40 60 アクチン骨格構造の変化をリアルタイムで観察 Time, t (sec) することが重要となると考えられる.また,骨 図 9 生骨組織への局所変形付与における蛍光輝 組織の変形が骨細胞に伝わるためには,骨細管 度比の時間変化 /骨細胞細胞突起間に何らかの刺激伝達物質の 存在が予想される.よって,骨細胞の焦点接着班の変化も同時に観察する必要があると考えられる. 4. おわりに 本研究では,骨細胞ネットワークが骨組織に付与された力学刺激を感知するメカニズムを解明する第 一歩として,単一の単離骨細胞および生骨組織中の骨細胞の局所変形に対するカルシウム応答特性に着 目した.その結果,単一の単離骨細胞は,細胞体ではなく細胞突起に対する局所変形に対してカルシウ ム応答を発現させることが示唆された.また,生骨組織に対するマクロな局所変形に対するミクロな骨 細胞のカルシウム応答を初めて確認できた.今後,骨細胞の力学刺激感知メカニズムが解明されれば, 骨組織・細胞の適応的リモデリング・力学刺激センシングをひとつのシステムと捉え,将来的にスマー ト複合材料の創製原理を確立できると期待される. 本研究は,京都大学大学院工学研究科機械理工学専攻修士課程学生の伊藤慎一君および神戸大学大学 院総合人間科学研究科人間環境科学専攻修士課程学生の青沼有紀さんの協力を得て進められた.ここに 記して,謝意を表す. 参考文献 [1] Knothe Tate, M.L., Adamson, J.R., Tami, A.E. and Bauer, T.W., The osteocyte, The International Journal of Biochemistry and Cell Biology, 36, (2004), 1-8. [2] Kamioka, H., Honjo, T. and Takano-Yamamoto, T., A three-dimensional distribution of osteocyte processes revealed by the combination of confocal laser scanning microscopy and differential interference contrast microscopy, Bone, 28-2, (2001), 145-149. [3] Cowin, S.C., Moss-Salentijn, L. and Moss, M.L., Candidates for the mechanosensory system in bone, Journal of Biomechanical Engineering, 113-2, (1991), 191-197. [4] Pead, M.J., Suswilo, R., Skerry, T.M., Vedi, S. and Lanyon, L.E., Increased 3H-uridine levels in osteocytes following a single short period of dynamic loading in vivo, Calcified Tissue International, 43 (1988), 92-96. [5] Van der Plas, A. and Nijweide, P.J., Isolation and purification of osteocytes, Journal of Bone and Mineral Research, 7, (1992), 389-396. [6] Kato, Y., Windle, J.J., Koop, B.A., Mundy, G.R. and Bonewald, L.F., Establishment of an osteocyte-like cell line, MLO-Y4, Journal of Bone and Mineral Research, 12, (1997), 2014-2023. [7] Klein-Nulend, J., Van der Plas, A., Semeins, C.M., Ajubi, N.F., Frangos, J.A., Nijweide, P.J. and Burger, E.H., Sensitivity of osteocytes to biochemical stress in vitro, FASEB Journal, 9, (1995), 441-445. [8] Kamioka, H., Miki, Y., Sumitani, K., Tagami, K., Terai, K., Hosoi, K. and Kawata, T., Extracellular calcium causes the release of calcium from intracellular stores in chick osteocytes, Biochemical and Biophysical Research Communications, 212, (1995), 692-696. [9] Miyauchi, A., Notoya, K., Mikuni-Takagaki, Y., Goto, M., Miki, Y., Takano-Yamamoto, T., Jinnai, K., Takahashi, K., Kumagawa, M., Chihara, K. and Fujita, T., Parathyroid hormone-activated volume-sensitive calcium influx pathways in mechanically loaded osteocytes, Journal of Biological Chemistry, 275, (2000), 3335-3342. [10] Bacabac, R.G., Smit, T.H., Mullender, M.G., Dijcks, S.J., Van Loon, J.J.W.A. and Klein-Nulend, J., Nitric oxide production by bone cells is fluid shear stress rate dependent, Biochemical and Biophysical Research Communications, 315, (2004), 823-829. [11] Bakker, A., Klein-Nulend, J. and Burger, E., Shear stress inhibits while disuse promotes osteocyte apoptosis, Biochemical and Biophysical Research Communications, 320, (2004), 1163-1168. [12] Thi, M.M., Kojima, T., Cowin, S.C., Weinbaum, S. and Spray, D.C., Fluid shear stress remodels expression and function of junction proteins in cultures bone cells, American Journal of Physiology-Cell Physiology, 284, (2002), C389-C403. [13] Kamioka, H., Sugawara, Y., Honjo, T., Yamashiro, T. and Takano-Yamamoto, T., Terminal differentiation of osteoblasts to osteocytes is accompanied by dramatic changes in the distribution of actin-binding proteins, Journal of Bone and Mineral Research, 19-3, (2004), 471-478. [14] Adachi, T., Murai, T., Hoshiai, S. and Tomita, Y., Effect of actin filament on deformation-induced Ca2+ response in osteoblast-like cells, JSME International Journal, C, 44-4, (2001), 914-919. [15] 田中基嗣,青沼有紀,安達泰治,上岡寛,山本照子,北條正樹,力学刺激に対する単離骨細胞の応 答観察,日本機械学会 2005 年度年次大会講演論文集, 05-1, 6 (2005), 119-120. [16] 田中基嗣,青沼有紀,安達泰治,上岡寛,山本照子,北條正樹,メカノセンサーとしての骨細胞に おける力学刺激応答の in vitro 観察,日本複合材料学会第 30 回複合材料シンポジウム講演要旨集, (2005), 255-256. [17] 青沼有紀,田中基嗣,安達泰治,上岡寛,山本照子,矢野澄雄,北條正樹,局所変形負荷に対する 単離骨細胞のカルシウム応答,日本機械学会第 16 回バイオフロンティア講演会講演論文集, 05-53, (2005), 7-8. [18] Tanaka, M., Aonuma, Y., Adachi, T., Kamioka, H., Takano-Yamamoto, T. and Hojo, M., Observation of calcium response of isolated osteocyte to localized mechanical stimulus, International Federation for Medical and Biological Engineering Proceedings, 12 (The 12th International Conference on Biomedical Engineering), (2005), CD-ROM. [19] 伊藤慎一,田中基嗣,安達泰治,上岡寛,山本照子,ニワトリ胚頭蓋冠培養組織内における骨細胞 の力学応答観察,日本機械学会第 18 回バイオエンジニアリング講演会講演論文集, 05-66, (2006), 19-20. [20] 伊藤慎一,田中基嗣,安達泰治,上岡寛,山本照子,北條正樹,生骨組織への力学刺激に対する骨 細胞カルシウム応答のリアルタイム解析,Preprints of Fuji 21c Forum, (2006), 47. [21] Adachi, T., Sato, K. and Tomita, Y., Directional dependence of osteoblastic calcium response to mechanical stimuli, Biomechanics and Modeling in Mechanobiology, 2-2, (2003), 73-82. [22] Grynkiewicz, Poenie, M. and Tsien, R.Y., A new generation of Ca2+ indicators with greatly improved fluorescence properties, The Journal of Biological Chemistry, 260-6, (1985), 3440-3450. MgO 複合材薄膜の作成と評価 工学研究科機械理工学専攻 田中 義和 Abstract: A protective layer, which supplies secondary electrons to discharge space, is one of key devices for a energy conservation of the AC type Plasma Display Panels. We investigated the properties of Eu2O3 doped MgO composite thin films deposited by an electron beam deposition method for the protective layer. The intended evaporation source of the composites were prepared and composition that calculated with following equation; [Eu/(Mg+Eu)]. The range of [Eu/(Mg+Eu)] ratio of the evaporation source were from 0 to 0.02. Composition, surface characteristics, crystallinity and ion induced secondary electron emissivility (γi) of the films were measured using Rutherford backscattering spectroscopy (RBS), Atomic Force Microscopy (AFM), X-ray diffractometer (XRD), and a self-developed analyzer. The experimental results shows that the highest value of γ was observed between 0.002 and 0.005 of [Eu/(Mg+Eu)] ratio. Moreover in the MgO composite thin films, the γ increases with the increase in the Root Mean Square roughness (RMS-roughness) Keyword: plasma display panel, protective layer, ion induced secondary electron emissivity 1. はじめに 昨年,フラットパネルディスプレイ型テレビ(FPDTV)の国内販売台数はブラウン管型テレビの販売台 数を抜きシェア50%以上となった.応答速度が速く且つ視野角が広いこと,また大型化が容易という特徴 からプラズマディスプレイテレビ(PDPTV)は、大型化テレビの本命と位置付けられて来た.しかし,昨 年後半より始まった37”超の画面を持つ大型の液晶テレビ(LCDTV)の本格的な市場参入により、2005 年第4四半期にはシェアを若干落としている.他のFPDTVと競争のためには、PDPTVは弱点とされている 高精細化と発光効率改善という二つの課題を克服しなければならない.このうち画面の高精細化は放電セ ルの小型化を意味し,小型化は発光効率の低下をもたらすことより発光効率改善の問題にも行き当たる. PDPではプラズマ放電を利用し光を得るが,放電の発生及び維持のためには放電空間に電子が多く存在す ることが望ましい.保護膜はイオン衝撃により放電空間に二次電子を放出する役割を持ち,PDPの発光効 率改善の鍵となる.この保護膜には酸化マグネシウム(MgO)膜が使われており,イオン衝撃二次電子放出 比γiに関する研究結果より,γiは結晶方向や表面形状に依存するとされている1,2).しかし絶縁体の二次電子 放出の機構には不明な点が多く,駆動回路・セル構造などに比べ改善が遅れている.1990年代後半より MgOと金属酸化物,例えばチタニア(TiO2)など,の混合物を用いてγiの向上を試みた報告がなされている2). これらの研究では蒸着材と薄膜で添加物の濃度差が激しく投入電力もMgO単体に比べると大きくなり実 際の生産には使用し難いという欠点がある.この問題点を解消するため、添加物濃度を2%までとしたMgO 混合物の研究を行っており、今回は添加物に酸化ユウロピウム(Eu2O3)を用いた蒸着材料を試作し実験を行 った.作成した薄膜はラザフォード後方散乱法(RBS),X線回折法,原子間力顕微鏡(AFM),二次電子放出 比測定装置で,組成,結晶方位,表面形状,γiをそれぞれ測定した. 2. 薄膜試料の調整及び評価 2.1 薄膜試料の調整 原材料は純度 99.9%以上の粉末を使用した.MgO と Eu2O3 の混合は,計算式 Eu/(Mg+Eu)を用いて金属元素中に Eu が 0.001,0.005,0.010,0.020 存在するように配合し,ボールミルを 用いて混合・撹拌を行った.出来上がった混合物は造粒後,回転式打錠機にて成型を行い,大気雰囲気下 1400℃以上で焼結し直径約 5mm,厚さ約 2mm の円盤状の蒸着材料を得た.焼結物の密度は計算より求め たバルク密度に対して 90%以上であった.保護膜の作成は一般に電子ビーム蒸着法を用いて行われ,製膜 時には基板過熱及び蒸着チャンバー内への酸素導入がなされる.このため本研究でも電子ビーム蒸着法を 用い,製膜時に酸素導入を行った.基板加熱はためハロゲンランプを用いチャンバー内を加熱した.製膜 時の条件を表 1 に示す.基板には Si(100)及びカーボン基板を用い,Si 基板は蒸着前にフッ酸及びアセトン を用いて超音波洗浄を行い Si 基板表面の SiO2 及び油脂成分の除去を行った.また製膜した Si 基板は大気 雰囲気下 400℃2 時間のアニール処理を行なった. 2.2.評価 薄膜の物性評価のため,組成分析,表面形状,結晶配向及び γ 測定の 4 項目の測定を行った. まず,ラザフォード後方散乱法 (RBS;Tandem accelerator Model14117,SEIKO Instrument Inc.;イオン種:He2+, ビームエネルギー2.0MeV) を用い薄膜の組成分析を試みた.薄膜の表面形状は,原子間力顕微鏡 (AFM;Nano scope III,Digital Instruments;Tapping Mode )を使って測定し、得られたデータより平均二乗粗さ (RMS 粗さ)を求め表面形状の指数とした.結晶方位は X 線回折法[XRD;Rigaku RINT-Ultima III;X 線源: を用い調べた. 次にイオン衝撃二次電子放出比γi を二次電 Cukα線(モノクロメータ使用)θ-2θ法] 子放出比測定装置により測定した.γi の測定条件を表2に示す.イオン衝撃型の二次電子放出の機構は入 射イオンの運動エネルギーの影響を受けないポテンシャル型とイオンの運動エネルギーに影響されるカイ ネティック型に分類される.MgO の場合両者の閾値は 200eV とされる 3). PDP の放電セル内でのイオン の運動エネルギーは 50eV 以下 4)とも言われ,ポテンシャル型である.一方,本研究で用いたイオンの運動 エネルギーは 500eV であり,このイオンエネルギー域では MgO 膜より放出される二次電子の数はエネル ギーの強弱に影響されると考えられる.二次電子放出比γi はイオンによっても異なり、P. Riccardi らは Ar イオンビームを用いた場合, 照射イオンの運動エネルギーによるγi の変化量は, He イオンや Ne イオンを 照射ビームとして用いた場合に比べて非常に小さいことを報告している 5).そこで本研究では Ar+イオンを 測定に用いた.また酸化物膜の表面に起こる帯電の影響を減少させるため,薄膜に照射する Ar+イオンの 電流量を昨年より 1 桁少ない数百nA 台とした. 表 2 イオン衝撃二次電子放出比γI 測定条件 表 1 薄膜の蒸着条件 Base pressure Working pressure Partial pressure of O2 Temperature in the Chamber 2x10-3 Pa 3.2x10-2 Pa 3.0x10-2Pa 250 degree C Base pressure Working pressure Ion species Accelerate Voltage 2x10-3 Pa 7.5x10-3 Pa Ar+ 500eV 3.測定結果 蒸着材に用いた,MgO は昇華型であり Eu2O3 は溶融型である.電子ビーム蒸着法では物質を蒸発させ 基板に堆積させるため、物混合材料の蒸気圧が大きく異なると,蒸着材料中と膜中でその割合が異なる恐 れがある.そこでラザフォード後方散乱法(RBS)にて膜の組成調べ結果を表 3 に示した.Eu2O3 は MgO よ り真空中での蒸気圧が低く,蒸着を行うと膜中の Eu の割合が蒸着材料中よりも少なくなると予想された. しかしながら計算式 Eu/(Mg+Eu)で求めた Eu の割合は 0.01 までは材料と薄膜との間に比較的差が無いこと が判明した.しかし材料中の Eu の割合が 0.02 の場合には膜中の Eu の割合が 0.013 と大きく異なった.ま た膜組成中に占める酸素の割合は平均で 52.3%と理論値より 2%ほど割合が高かった. 次に AFM による観察結果を述べる.プローブには Nanoworld 社の NCH-10T シリコン単結晶プローブ を用いた.接触式の測定では針の先端が磨耗すると得られる画像イメージがいびつになり RMS 粗さが小 さくなる傾向がある.そこで通常のコンタクト式の測定に比べ針先の磨耗が少ないタッピングモードで測 定を行った.今回,鮮明な MgO 膜の RMS 粗さを調べてところ 3nm 以上であった.その為標準試験片と して MgO 膜を用い,サンプルの測定の前後に MgO 膜を測定し RMS 粗さが 3nm 以上であることを確認し ながら測定を行った.結果を図-1 に示す.Eu 濃度が高くなるにしたがって RMS 粗さの値が小さくなる傾 向を示している.しかし MgO 膜の RMS 粗さは一番粗いにもかかわらずγi が低くこのことは単に RMS 粗 さがγi を決定しているのではないことを示している. XRD の測定結果を図 2 に示す.γi 測定用の薄膜と同じ膜厚のサンプルを用いて測定を行った.何れの膜 でも(111),(200),(220)のピークが確認されたが,膜が薄いためピークの強度が弱くそれ以上の分析は出来 なかった.また膜中の Eu の割合が一番高い薄膜でも Eu2O3 のピークは確認出来なかった. 次に γi の測定結果について述べる. γi は薄膜中の金属原子に対し Eu の割合が 0.002 の場合に最大値を 観測した。 実験で得られた Zn の割合とγi より二次曲線の回帰式を求めると 0.002~0.007 の間で極大となっ た.その際のγi の値は MgO 膜に比べが約 10%高くなった.膜中に含まれる Eu の割合がそれ以上高くなる .極大値を取る添加剤の とγi は低下し Eu の割合が 0.013 の場合は,無添加の試料と略等しくなった(図 3) 割合は Ni または Zn を添加した場合と似ていることから,MgO 中へ金属酸化物の微量添加をする場合, 膜中の全金属元素中に添加物の金属の占める割合が0.002~0.005 にすることがγi の向上には望ましいと考え られる 6,7). 表 3 膜中の元素割合 Concentration of elements (at.%) Eu/(Mg+Eu) ratio of evaporation source Mg(%) Eu(%) O(%) Eu/(Mg+Eu) ratio of thin films 0.000 0.001 0.005 0.010 0.020 48.8 48.0 46.0 47.6 46.7 0.1 0.3 0.5 0.6 51.2 51.9 53.7 51.9 52.7 0.000 0.002 0.006 0.010 0.013 3.400 3.430 3.200 3.254 3.139 3.000 2.800 2.600 2.400 2.312 2.200 Eu割合 (200) 0.005 0.010 0.015 0.010 0.006 0.002 0.000 34 44 54 64 2θ (degree) Cation Ratio [Eu/(Mg+Eu)] 図 1 膜中の金属元素に含まれる Eu の割 合と膜の RMS 粗さの関係. (220) 0.013 2.189 2.000 0.000 (111) Intensity (arb.units ) RMS-Roughness (nm) 3.600 図2 MgO 混合物薄膜の XRD パターン. 4.まとめ 1. 2. 3. 4. 5. MgO と Eu2O3 の混合物の製膜を行いγi の測定を行った. 計算式 Eu/(Mg+Eu)で求めた Eu の割合が 0.002 のとき最も高いγi が確認された.また測定結果は 0.002~0.005 付近に極大値を持つ二次曲線に良く一致した.添加物の割合が 0.002 近辺でγi が高く なる傾向は本研究室で行われた他の研究でも見られ,この付近に最適な添加割合があると推測され る. 薄膜の組成分析の結果,蒸着材に含まれる金属元素の内 Eu の割合に関して,0.01 までは薄膜は蒸 着材の組成を反映していることが確認出来た.しかし Eu の割合が 0.02 のときには膜中の Eu の割 合が 0.013 と異なった.このことより Eu を添加する割合は 0.01 までが良いと考えられる. RMS 粗さは膜中の Eu の量が増えるにしたがって小さくなりγi も小さくなる傾向を示した.しかし MgO 膜の RMS 粗さが一番粗いにもかかわらずγi は低いことから単に RMS 粗さがγi を決定してい るのではないことを示している. 試験に用いた薄膜は(111),(220),(200)の結晶方位を持つことを確認した. Secondary electron emission (γ) 0.130 0.128 0.125 0.126 0.120 0.115 0.116 0.113 0.110 0.105 0.108 0.100 0 0.002 0.004 0.006 0.008 0.01 0.012 0.014 Cation ratio [Eu/(Mg+Eu)] 図 3 金属原子中の Eu の割合と二次電子放出比の関係 参考文献 [1] Y.H.Cheng, H. Kupfer, F. Richter, H. Giegenack, and W. Hoyer, “Structure and secondary electron emission properties of MgO films deposited by pulsed mid-frequency magnetron sputtering”, Journal of Applied Physics, 93 (2003), 1422-1427 [2] Y. Kim, R. Kim, H.J. Kim, H. Jeon, and J.-W. Park, “Effect of TiO2 Addition on Secondary Electron Emission and Discharge Properties of MgO Protective Layer”, Materials Research Society, Symp. Proc. 621, (2000), Q5.6.1-6 [3] 内池平樹・内田龍男, “フラットパネルディスプレイ大事典”, (2001),504,工業調査会. [4] K.S. Moon, J. Lee, and K.-W. Whang, “Electron ejection from MgO thin films by low energy noble gas ions”, Journal of Applied Physics, 86,7,(1999),4049 [5] P. Riccardi, M.Ishimoto, and R. A. Baragiola, “Ion-induced emission from MgO by exciton decay into vacuum”Surface Science Letter,571,(2004),L305 [6] 村上茂人, 井手亜利,他 “MgO-NiO 複合材料薄膜の二次電子放出比の測定”,日本機械学会関西支部 第 80 期定期総会講演会予稿集(2005), 9-23-24 [7] 田中義和, “PDP 保護膜材料開発に関する研究”, 21 世紀 COE プログラム「動的機能機械プログラム」 平成 16 年度年次活動報告書, (2005), 179-182 静電型 MEMS デバイスを用いたナノ材料の電気機械特性評価 工学研究科マイクロエンジニアリング専攻 土屋 智由 Abstract: This report describes the research on a new on-chip tensile testing device for nano-scale materials such as single walled carbon nanotube. This device consists of a comb drive electrostatic actuator for tensile force application, two sets of parallel plate capacitance for force and elongation measurement, and a folded-beam spring to convert the tensile force into the displacement. The device was fabricated from silicon on insulator (SOI) wafers using surface micromachining processes. The key for high resolution and accuracy of the force measurement is the precise measurement of the folded-beam spring structure. The spring constant was measured using a tensile tester with electrostatic grip system. The force resolution of the tester is less than 0.2 µN. The spring constant of two type beam structure was 0.18 N/m and 0.57 N/m. Key words: mechanical properties, electrostatic transducer, MEMS, carbon nanotube. 1. はじめに カーボンナノチューブを代表としたナノスケールの材料(以下ナノ材料と書く)はサイズ効果や表面 効果,量子効果などの発現によって特異な物性を示すことが期待されている.機械材料としてもナノ材 料は高い弾性率,強度などを持つことが数値計算によって予想され,複合材料や微小な機械構造体 (MEMS; Micro Electro Mechanical Systems)などへの応用が検討されており,ナノ材料を含む複合材料 の機械的特性の測定からその優れた性質が明らかになりつつある.しかしながら,ナノ材料単体として の機械的特性を直接測定することは寸法,測定すべき物理量の微小さゆえに非常に困難である. 本研究ではナノ材料単体の機械的特性の測定のため,材料試験において基本的かつもっとも単純な測 定法である単軸引張試験を実現することを目指している.たとえば単層カーボンナノチューブ (SWCNT)はヤング率が 1 TPa,強度が数十 GPa という値が報告されており[1],ここから長さ 1µm の SWCNT の剛性は 1 N/m,破断荷重が数十 nN,破断伸びが数十 nm ということが予想され,既存の報告 されたマイクロ材料の引張試験法では測定不可能である.そこで本研究では数 nN,数 nm の荷重,変位 の計測が可能な静電容量による変位検出に着目した.微小な静電容量の検出分解能は実用化が進む加速 度センサなどにおいて実証されており,0.1nN,0.1nm の分解能が十分可能である.本報告では静電容量 型引張試験デバイスの設計と試作について述べた後に,荷重を変位に変換する可撓体として用いた折り 返し梁構造の剛性を機械的測定法により較正した結果について報告する. 2. 静電容量型引張試験デバイス 2.1 原理 本研究で提案する引張試験デバイスの構造を図1に示す.膜厚が数~20 µm程度の膜構造で表 面マイクロマシニングを用いて作製されることを意図している. SWCNTなどのナノ材料の試験片は下部 のギャップに固定され,櫛型静電アクチュエータによって試験片上端を上方に変位させることで引張荷 重が印加される.試験片に印加される荷重はデバイス中央にあるバネで変位に変換される.試験片の伸 びとバネの変位はそれぞれ下部,上部に設置された平行平板静電容量(キャパシタンス)で検出する. 可動部全体は4本の梁で支持される.これらの支持梁はデバイスの上部(荷重-変位変換用ばねよりも上) に取り付けられているのでその剛性は荷重検出には影響しない. 2.2 設計 SWCNT の測定を考慮して,デバイス各部の設計を行う.SWCNT のヤング率を 1TPa, 長さ を 1 µm とし,断面積が 10-19~10-18 m2 であるとすると試験片のバネ定数 kSAMPLE は 0.1~1 N/m と予想され る[1].また,強度が数十 GPa であることから試験片の破断ひずみが数%,すなわち伸び,破断荷重はそ れぞれ,数 10nm,1~数 10 nN 程度になる.測定においては伸び,荷重をそれぞれ破断時の 1%の分解能 で測定することを目標とする. ここでデバイスの特性を考える.まず,伸び検出部の変位 d1 は d1 = FSAMPLE k SAMPLE (1) で示される.試験片に加わる引張荷重 FSAMPLE と等しい力が荷重変位変換ばねに加わり,バネは変形す る.荷重検出部の変位 d2 は d 2 = FSAMPLE k + d1 = FSAMPLE (1 k + 1 / k SAMPLE ) (2) である.静電アクチュエータの駆動力 Fcomb は引張荷重と可動部支持用梁の復元力 FSB の和であるので, Fcomb = FSAMPLE + k SB d 2 = FSAMPLE (1 + k SB / k + k SB / k SAMPLE ) (3) である.式 2 より荷重検出部の変位は試験片の伸びとばねの変形の和であるのでばねの変形を精度よく 測定するためにはばね定数 k が試験片の剛性と同程度かそれより小さいことが望ましい.剛性の低いば ねを作製することは作製歩留まりを考えると困難であると考え,k = kSAMPLE とする.一方,支持梁の剛 性 kSB は作製歩留まりを考えると大きいほうがより望ましい,しかし櫛型静電アクチュエータの発生力 を考慮して最適値を検討する必要がある.以上の点を考慮してデバイスを設計した. 引張荷重印加: 櫛型静電アクチュエータ 可動部支持用梁 バネ定数: kSB 変位: d2 荷重検出: 平行平板静電容量 固定 ビーム長さ: L1 バネ定数: k 荷重-変位変換可撓体: 折り返しばね 伸び検出: 平行平板静電容量 固定 変位: d1 試験材料 バネ定数: kSAMPLE 図 1 デバイス構造 表 1 荷重-変位変換ばねの寸法 2.3 作製プロセス 作製プロセスが単純な SOI ウエ 0.3 0.5 0.7 1.0 ハを用いた単結晶シリコン構造により上記のデバイス ばね定数 (N/m) 構造を実現する.試験片である SWCNT の形成を考慮 2 幅 (µm) 283 229 214 190 して構造の厚みを 5 µm と設定して具体的な寸法を設 長さ (µm) 5 膜厚 (µm) 計した.結果,後述の剛性測定で用いるばね部につい ての詳細な寸法を表1に示した.作製プロセスはワイ ヤボンディングパッドのアルミ電極とデバイス構造の 2 枚のマスクを用い,アルミ電極を加工した後に ICP-RIE によるドライエッチングで単結晶シリコンにデバイス構造を作製する.その後 SOI ウエハの埋 め込み酸化膜を犠牲層としてバッファードフッ酸でエッチングし,構造体をリリースする. 3. 微小ばね構造の剛性測定 上述のナノ材料引張試験デバイスの測定分解能,精度を決定する重要な要素の一つとして引張荷重を変 位に変換する折り返しばね構造が挙げられる.このばねの剛性すなわちばね定数を正確に知ることで引 張試験デバイスの荷重測定の精度を高めることができる.しかしながら,このような膜厚数µm,長さ数 百µm のばねの剛性は共振周波数などの間接的な測定のみが行われており,直接引張荷重を加えてその 荷重と変位を測定することはほとんど行われていない.われわれは微小な薄膜材料を引張試験する装置 を開発し,膜厚数十から 100nm 程度の試験片を試験するための装置を製作,評価している[2].そこでこ こではこの装置を用いて上記のばね構造の荷重-伸び線図を 測定することを試みる. 3.1 測定法 ばね定数の測定には引張試験デバイスと同時 に作製したばね定数評価デバイスを用いた. (図 2)上記のば ね構造を 4 個並列に接続し,その一端を基板に固定し,他方 は引張試験装置の荷重軸へのチャッキングのため,矩形の平 板になりここは基板から犠牲層エッチングでリリースされる. この試験片の自由端は静電チャック法で固定する. 3.2 測定装置 試験装置の外観と構成図を図に示す.荷重測 定機構は熱重量測定器(島津製作所 TGA-51H)の荷重測定部 を用いた.これはサーボ型の天秤になっており,釣り合いを 保つように電磁コイルに電流を流し,この電流を荷重として 読み取る.荷重測定端(図中左端)を荷重軸としてここに静 電チャック機構を接続した.試験片チップは X, Y, Z に変位可 図 2 ばね定数評価デバイス 能な試料台に固定されていて Z 軸(引張軸)と X 軸(チャッ クと試験片の間隔)方向にはモータ駆動のステージを用いている.また,Z 軸ステージの変位はレーザ ー変位計を用いて測定した.本試験では被測定部(ばね)以外の剛性が被測定部よりも十分高いのでス テージの変位をばねの変形と同じとみなす. Force coil Weight measurement unit of Shimadzu TGA-51H analyzer Electrostatic grip Counter weight Specimen Z Force-feedback current Picomotor Displacement Photoelectric transducer To computer Control circuit Force X, Y stage Stage driver Displacement 図 3 静電チャック方式の微小荷重引張試験装置 3.3 結果と考察 試験時のばねの変形を約 5 秒おきに撮影した(図 4) .試料台ステージを変位させてい るので,図では固定部が下方に移動し,静電チャック(図上方)は変位しない.4 つのばねが等しく, また滑らかに変形していることがわかる. ① ② ③ ④ ⑤ ⑥ ⑦ ⑧ 500µm 図 4 引張試験時のばねの変形(5 秒おきに撮像) 250 250 200 200 150 150 荷重(µN) 荷重(µN) 図 5 に荷重-変位の測定結果を示す.荷重増加時,減少時のプロットはよく一致している.また,繰 り返しの測定による変動もない.これらのプロットの傾きから得られた 1 個のばねの剛性はそれぞれ, 0.18 N/m (設計値 0.3 N/m),0.57 N/m (設計値 0.7 N/m)であった.設計値との乖離の原因はばねの断面形状 にあると考えられる.断面の SEM 観察によると本来矩形であるべきばねの断面が基板側の下部が大き くえぐられた形状になっていた.この断面形状を考慮して剛性を計算したところそれぞれ 0.25 N/m,0.56 N/m となり,これらは測定値に近い値であり,計測値は妥当であると考える.このように解析,FEM な どによる剛性の計算では十分予測できない微小なばねの剛性を直接測定できることを確認した. 100 50 50 0 0 -50 -20 100 0 20 40 60 変位(µm) 80 100 -50 -20 0 20 40 60 80 100 120 変位(µm) 図 5 ばねの荷重-変位線図(4 本のばねの合計) 4. まとめ カーボンナノチューブなどのナノスケール材料の機械的特性を測定するための静電型 MEMS 構造を用 いた引張試験デバイスを提案し,デバイス構造を設計,試作を行った.また,サーボ型天秤を用いた静 電チャック型引張試験装置を用いて,このデバイスの荷重測定の精度を決定するばね構造の剛性を機械 的に直接測定した.今後は,試験片を引張試験デバイスに集積化する手法を検討し[3],ナノ材料の引張 試験を実現する.さらに MEMS デバイスの特徴を生かし,電気的特性を同時に計測するデバイスの実 現を目指す. 参考文献 [1] V. N. Popov,Carbon nanotubes: properties and application,Mat. Sci. Engineering R43 (2004), 60-102. [2] T. Tsuchiya, M. Shikida, and K. Sato, Tensile testing system for sub-micrometer thick films, Sensors and Actuators A 97-98 (2001), 214-228. [3] C. Stampfer, A. Jungen and C. Hierold, Single walled carbon nanotubes as active elements in nano bridge based NEMS, Proc. 5th IEEE Conference on Nanotechnology, Nagoya, Japan July 11-15, 2005, TU-P8-3. 磁場中粒子挙動解析および構造制御に関する研究 工学研究科マイクロエンジニアリング専攻 津守 不二夫 Abstract : It is important to control micro structure of the powder particles, which could be used to make new functional materials. One of effective methods to arrange particles is application of magnetic field. The magnetic particles form chain or columnar structures with the applied magnetic field. To design the structure, a new computational system was developed to simulate this magnetic powder behavior. Discrete Element Method (DEM) was employed and coupled with a magnetic FEM to estimate magnetic field. Using the present system, growth of column powder structures is simulated and compared with experimental ones. It is found that the width of the columnar structures changes with magnitude of applied magnetic field in both simulated and experimental results. Key words: Discrete Element Method, DEM, powder material, magnetic field, micro actuator. 1. はじめに 磁性粉末は磁性粒子の集合体であり,外部磁場の印加により柱状構造等の微細な構造が生成す ることが知られている.このような微細構造は希土類焼結磁石材料の性能や磁性流体の特性に影 響を与える.さらには,このような鎖状構造を利用し MEMS 等への微細構造の作製の応用が期待 されている.粉末挙動を把握することは非常に重要であり,磁場中の粉末の挙動を実験的・解析 的に把握する試みが近年行われている.粉末粒子の挙動を解析する手法として,個別要素法 (Discrete Element Method; DEM)の利用が計算機の発達とともに広がっている[1].これは分 子動力学法をもとにした解析手法であり,粉末粒子のひとつひとつの動きをすべて評価する方法 である.粒子はお互いに衝突することにより反発し,また接する場合には摩擦の影響を考慮する 等,ミクロな粉末同士の相互作用を明確に解析に導入することができる.磁場中での粒子構造の 把握には粒子のひとつひとつを直接追跡する必要があるため,本研究ではこの DEM を利用するこ ととする.DEM での粉体挙動の解析は,圧縮成形や粉末流動のシミュレーションに応用されてお り報告例も多い.特に磁場の印加に起因する柱状・鎖状構造の発生のような現象は連続体力学的 な手法での解析の糸口はまったくつかめておらず DEM での解明が期待されており,実際に DEM に おいて鎖状構造の生成解析例が報告されている.しかし,従来の報告例には問題も多い.報告さ れている解析例では外部から与えられた磁場による電磁力が粒子を移動させているが,ここで問 題となるのは磁場の取り扱いである.解析内では簡単のため,磁場の影響に大胆な近似を組み込 まれている.特に粒子の磁気モーメントがその粒子の中心に集中するという仮定[2]が広く使用 されている.しかし実際には粉末内部の磁場の状態は粉末粒子が存在するだけで乱れが生じ,こ こで発生した乱れは粒子の移動にも影響するため,このような乱れを正しく評価することが重要 である.さらには一般に磁性材料の B-H 関 係は線形関係にないため単純な仮定のみで は問題が発生する可能性もある.この問題 Geometry に対処するために個々の粒子に関しても全 of particles てミクロレベルでメッシュ分割した磁場 FEM 解析を DEM 解析に連成させたシステム Magnetic を完成させた[3-6].磁場 FEM 解析において force は非線形な磁気的材料特性も導入が容易で Magneto-FEM Distinct Element Method あり,磁場の影響のより厳密な取り扱い が可能である. 図 1 FEM-DEM 連成モデル. このような解析システムとリンクする ことにより粉末材料の磁場による影響を明らかにし,また構造の設計へと応用が可能となる.磁 場の利用は直接個々の粒子をマニピュレートすることなく,全体構造にマクロに外部磁場を与え ることにより粉末粒子レベルの微細な構造を生成できるため,生産性の高い有効な手法である. これにより機能性材料の生産も可能となる.例としては高性能焼結磁石や磁性流体がこれにあた る.この他にミクロ構造自体を微細部品として取り扱うこともできる.本研究では外部磁場によ り磁性粒子構造を制御し,弾性樹脂内部に粒子構造を固定したマイクロアクチュエータを作製す ることを試みる. 2. 解析方法 有限要素法による磁場解析と DEM の連成解析システムを作成する.これを利用して外部磁場が かかった場合の粒子挙動や,壁面を移動させた場合の構造変化について解析を行う. 2.1 連成方法 解析においては磁場中成形プロセスを微小なタイムステップに分け,そのタイム ステップごとの粒子の加速度・速度・位置を追跡することにより各粒子の状態の時間履歴を得る. まず,各タイムステップにおいて有限要素法による磁場解析を行い,個々の粒子に作用する電磁 的な力を求める.この力を DEM に取り込むことにより粒子個別要素法と有限要素磁場解析法の連 成を行う.DEM での粒子位置を導出する際には発散を防ぐため有限要素解析の 1 ステップに対し, DEM では 1000 ステップと十分に小さい 1000 倍の時間分割を行った. 2.2 有限要素磁場解析 本研究では,磁場解析ソフト FEMM を用いて,有限要素磁場解析を行う. これにより,解析領域各要素での磁場の向き,磁束密度,各粒子にかかる力及びポテンシャルエ ネルギーを算出する.FEMM における磁場解析では A をベクトルポテンシャルとして,以下に示す ラプラスの方程式を支配方程式とする. ∇2 A = ∂2 A ∂2 A + =0 ∂x 2 ∂y 2 (1) ある要素 p におけるベクトルポテンシャルを, A = α 0 + α1 x + α 2 y と近似する. る. (2) α 0 , α 1 , α 2 は係数である.このときの残差を r として,式(1)は以下のようにな ∇2 A = ∂2 A ∂2 A + =r ∂x 2 ∂y 2 (3) この誤差 r が最小となるような各要素におけるベクトルポテンシャル A を求め,ここから磁束 B を求め,各要素の持つエネルギーが導出される.さらに Maxwell の応力テンソルから各要素に働 く力を求めることができる. 3. 粒子配列 PDMS 樹脂中に粒子を分散させ,磁場中で 樹脂を硬化することにより粒子構造を樹脂中 に固定する.用いた粒子は 64~73 μm に分級 された鉄粉末であり,図 2 に示すように樹脂 に分散させた後スライドガラス,カバーガラ スにはさみこみ 100 μm の厚さとなるように 硬化した.粒子は樹脂との重量比で 0.1~0.3 PDMS & particles Defoaming Slide glass Apply magnetic field 図 2 サンプル作製法. Cover glass 1 mm -2 2.08×10 T 1 mm -2 5.20×10 T Applied magnetic field 2.08×10-2 T 5.20×10-2 T 図 3 磁場中粒子挙動(左:解析結果,右:光学顕微鏡観察結果) で混合した.これをコイルの中に設置し,コイルにより発生させた 2.08×10-2~5.20×10-2 T の磁場 中で硬化した.磁場はスライドガラスの面内方向とした.得られた樹脂中の構造は光学顕微鏡に より観察する.これらの実験結果は前節の手法による2次元解析結果と比較する. Distance(μm) 結果 図 3 に解析結果および実験結果を示す.両者ともに外部磁場の方向に並行な鎖状・柱状の 構造が生成している.柱状の構造は外部磁場の大 600 きさにより構造が異なっていることがわかる.つ まり,外部磁場が増大するにつれてより太い構造 となり,同時にそれぞれの構造間隔も広がってい 400 る.図 4 は外部磁場強度と柱状構造間の平均距離 の関係をプロットしたものである.磁場の増大と ともに距離が増大している関係が見てとれる.間 0.1 experiment 隔は粒子密度の高いものほど小さくなっており, 200 0.2 experiment 柱状構造の本数自体が増加することがわかる.こ 0.3 experiment のように柱状構造の太さが増大することは,粒子 0.2 simulation(water) の接触部位の磁束の飽和によるものであることは 0 0.02 0.04 0.06 既報のとおりである[3]. Flux density(Tesla) 前節に示したように外部磁場を利用すること により樹脂中に磁性粒子を配置した構造を容易 に作製することが可能である.このような磁性 粒子が肯定された構造に磁場を印加することに より粒子間に力を発生させ,シートを変形させ ることができる.弾性シート状にインクジェッ トにより規則配置した磁性粒子ドットを利用し シートを変形させることも可能である[7].しか し,この方法では粒子配置に自由度はあるもの の大量生産にはまったく向かない.本節ではこ のような短所を克服するために,柱状構造を有 する樹脂構造がそのままアクチュエータへと応 用可能であることを示す.図 5 はアクチュエー タの作製方法および原理である.3 節のサンプ 図 4 外部磁場強度と粒子による柱状構造間 距離との関係 PDMS+particles 2nd PDMS layer Magnetic field 4.アクチュエータへの応用 図 5 内部構造を含む弾性樹脂シートによる マイクロアクチュエータ Displacement(μm) ルと同等の,粒子を含む樹脂の上にもう一層粒子 を含まないシートを配置する.このようにして作 sample 1 100 製された2層のシートは磁場による伸縮挙動の違 sample 2 sample 3 いから,外部磁場により反り変形を示す. sample 4 シートの作製条件としては鉄粒子を重量比で sample 5 -2 0.1,外部磁場 2.08×10 T をかけ厚さ 100 μm のシー 50 トを作製し,30 μm の厚さの PDMS 樹脂シートを 積層させる.この積層シートを 1.0 × 5.0 mm に切 り出したものを図 5 の右に示すように上端を固定 し,コイル中に設置した.そして外部磁場を変化 させながら下端位置の変化をレーザー変位形によ 0 0.02 0.04 り測定した.図 6 に作製した 5 個のサンプルの測 Flux density(T) 定結果を示す.下端の変位量は磁場強度に応じて 図 6 マイクロアクチュエータの変形. 変化している.ただし同条件で作製したサンプル であるにも関わらず個体差の大きいこともわかる. これは粒子サイズと粒径分布のばらつきによるも のが原因であり,今後微細な粒子や単一粒径の粉末材料を利用することに改善が可能であると思 われる. 5. まとめ 外部磁場により磁性粒子構造を樹脂中に制御して配列させることが可能となった.このような 構造の生成について実験のみならず解析によりその構造を予測できることを示した.また,樹脂 中に磁性粒子構造を配置することによりアクチュエータを作製した.従来よりも容易にサイズダ ウンが可能であり,マイクロアクチュエータの作製方法として有望と思われる.内部の粒子構造 を変化させることにより,一様な磁場においてもさまざまな変形挙動を示すアクチュエータも設 計可能であり,今後磁場によるさらに複雑な構造制御が期待される. 参考文献 [1] Cundall, P. A and Strack, O.D.L , “A discrete numerical model for granular assemblies”, Geotechnique, 29 (1979) 47-65. [2] R.S. Paranjpe, “Stability of chains of permeable spherical beads in an applied magnetic field”, J. Appl. Phys. 60-1(1986) 418-422. [3] 津守不二夫, 平田正道, 島進, “FEM-DEM 連成モデルによる磁場中粒子挙動解析”, 粉体および粉 末冶金, 52-3(2005), pp194-198. [4] 津守不二夫,栗原史和,小寺秀俊,島進, “磁場中粉体成形における粉末流動のキャビティ形 状による影響”, 粉体および粉末冶金, 52-6(2005), pp458-463. [5] Fujio Tsumori, Masamichi Hirata and Susumu Shima, “Column Structure Growth Simulation of Magnetic Particles by Distinct Element Method Coupled with Magneto-FEM”, Proc. Powder Metallurgy & Particulate Materials, vol. 1, (2005), pp 74-84. [6] Fujio Tsumori and Susumu Shima, “Simulation of Column Structure Growth of Magnetic Powder in Applied Magnetic Field by Coupled FEM-DEM Modeling”, Proc. Intelligent Processing and Manufacturing of Materials (2005). [7] Fujio Tsumori, Hiroyuki Ogawa and Susumu Shima, “On-demand Droplet Spotter for Formation of Micropattern Using Electrostatic Force”, roc. Powder Metallurgy & Particulate Materials, vol. 9, (2005), pp 96-107. a-SiO2 のモデル化と外部電磁場に対する電子状態の応答 についての理論研究 工学研究科マイクロエンジニアリング専攻 土井 謙太郎 Abstract: A modeling procedure of amorphous structures has been studied using molecular dynamics (MD) simulation and first-principle calculation. In the MD simulation, temperature and pressure were controlled and then lattice deformations were considered. Meta-stable structures and electronic structures of amorphous silicon dioxide (a-SiO2) were investigated. The electronic structures showed some differences between the amorphous structures and pure crystal. For example, in amorphous phases, the band edges of SiO2 were disturbed because of changes of the bond lengths and angles between Si and O atoms and generation of the dangling bonds. In this work, behavior of oxygen vacancies was also simulated and the electronic structures caused by the vacancies were investigated. These electronic structures were discussed by density of states (DOS) and differences between the pure crystal and amorphous phases were revealed. Some interesting structures were observed by our MD simulation and their chemical-bonding properties were clarified by quantum energy density. Key words: a-SiO2, Molecular dynamics、Electronic structure, First-principle calculation, Electromagnetic field 1. はじめに 近年、半導体の分野において、ロードマップに従ったデバイスの急速な微細化により,理論的な考察が不十 分なままものづくりが先行している感がある.そこには,ナノオーダのデバイス中において,これまでの古典的 な理論体系の枠組みでは理解できない現象が多々現れることが挙げられる.これまでは,シリコンを基盤にし た半導体デバイスが主流であったが,ロードマップに従えば,従来の方法による微細化ではシリコン半導体は 限界に達しているといえる.特に,電界効果トランジスタ(MOSFET)の信頼性の要であるゲート絶縁膜は 2nm を 切る膜厚が要求される.しかし,このオーダの絶縁膜ではゲート電極とシリコン基板間のリーク電流の抑制が極 めて困難とされる.そこで,現在では絶縁膜の膜厚を厚くする代わりに誘電率を上げる,添加物によりリーク電 流を抑制する,さらにはトランジスタの構造そのものを従来型から変えてしまうなどの試みがなされている.また, 実験グループにおいても,新規材料の探索のために第一原理計算を用いることが盛んになされるようになって きた.これまでのあらゆるものについて試行錯誤するという方法に比べ,コストと時間の面で大幅な削減が期待 されるからである.第一原理計算は経験的なパラメータを用いることなく Shrödinger 方程式を解くことにより 系の電子状態が計算できる利点があるが,現在の計算機の能力では 102 オーダの原子数を扱うのが限界 である.この制限の中で,現実の構造や現象をモデル化し考察することは,まだ発展段階であり,さま ざまな分野で議論がなされているところである. 我われは,これまで第一原理に基づく電子状態計算により,シリコン酸化膜について電子状態の考察 を行ってきた[1,2].本研究においては,MOSFET のゲート絶縁膜に注目し,その安定構造をモデル化し, さらにその電子状態を考察することにより,絶縁膜のリーク電流発生機構の解明を最終目標として研究 を進めている.その第一歩として,平成 15 年度は,固体中の伝導電子が占めるフェルミ面の計算を行う ことにより,結晶中の伝導に寄与する電子に注目した考察が可能となった.つまり,電流が流れると考 えられる電子構造が得られれば,その中の電流に寄与する電子のみに注目した議論が可能となる.引き 続き,今年度においてはシリコン酸化物(SiO2)のアモルファス構造(a-SiO2)をモデル化することにより, 現実の膜構造をモデル化した上で,その電子状態について議論した.具体的には,温度と圧力を制御し た分子動力学(MD)計算により,SiO2 の結晶格子を高温でアニールした後に急冷することで準安定なアモ ルファス構造をモデル化した.さらに,MD 計算で得られたアモルファス構造に対して第一原理計算を 行うことによりその電子状態を計算した.本研究においては,MD 計算には Vashishta らによって開発さ れた経験的ポテンシャルを用いた古典 MD 計算を行った.この手順により,絶縁膜の破壊の原因のひと つと考えられている酸素欠陥の挙動についてシミュレーションを行い,その電子状態について考察を行 った結果を示す.なお,本研究における計算には京都大学学術情報メディアセンターの大型計算機シス テム HPC2500 を利用した. 2. 計算方法 2.1 分子動力学計算 本研究では原子数,温度,および圧力を一定に保った MD 計算を行うために,ハミ ルトニアンを式(1)のように定義する. ( ) r& t t r& 1 wh ws s& 2 1 1 t& & H = ∑ mi ρ i h hρ i + ∑ ∑ φ (rij ) + 2 tr h h + Ω tr P + + 3 Nk BT ln s 2 i 2 i j ≠i 2s 3 2 s2 r ここで,i 番目の粒子の位置は ri = (1) r ρ i h と表され,h は基本並進ベクトルからなる 3×3 行列である.s は 温度を一定に保つための変数であり,Ω は体積,P は外部の圧力である.wh,ws はそれぞれ h,s の共役 量に現れる質量にあたる.本計算では適当なパラメータとして扱っている.さらに,N は全粒子数,kB は Boltzmann 定数,および T は外部温度である. 式(1)そのものは正準方程式を満たさないため,厳密 にはハミルトニアンと呼ぶことはできないが,適当な変換により正準方程式を導き出せる形に変形でき る.詳細については文献[3]に書かれている.また,式(1)から dH dt = 0 が確かめられることから,MD 計算をする上で H を運動に関する保存量とみることができる.式(1)は NVT 系を扱うために Nosé によ って考案された方法[3]と,NPH 系を扱うために Parrinello と Rahman によって考案された方法[4]から, 本研究において導き出したものである.式(1)から導出される運動方程式は以下のようになる. ( ) r s& r ρ&&i = − ρ& i − h t h s −1 r 1 ∂φ (rij ) ρ ij d t r& h h ρi − ∑ ∑ dt ∂rij rij i j ≠ i 2 mi ( ) rr 2 ⎛ ⎞ ( ) ∂ φ r r & r r 1 s s ij ij rij && = h& + ⎜ ∑ mi vi vi − ∑ ∑ − PΩ ⎟ h t h ⎟ s wh ⎜⎝ i 2 i j ≠i ∂rij rij ⎠ (2) ( ) &s& = −1 r r ⎞ s& s ⎛ + ⎜ ∑ mi ρ& it h t hρ& i − 3Nk BT ⎟ 2 s ws ⎝ i ⎠ (3) (4) r r ここで,i 番目の粒子の速度は vi = ρ& i h で定義している.φ(rij)は 2 粒子間のポテンシャルであり,本計算 ではVashishtaらによる3体間の効果を考慮したポテンシャル[5]を用いた. 外部の圧力と温度はそれぞれ, 式(3)と式(4)の散逸項より以下の条件を満足しなければならない. rr ∂φ (rij ) rij rij rr 1 1 ⎛⎜ P= ∑ mi vi vi − 2 ∑i ∑ ∂rij rij Ω ⎜⎝ i j ≠i T= 1 3Nk B rt r ∑ m ρ& h hρ& i t i ⎞ ⎟ ⎟ ⎠ (5) (6) i i 3 体間の相互作用を考慮した経験的ポテンシャルは以下のような関数で与えられる. φ (rij ) = ∑ φij(2 ) (rij ) + i< j r r ∑ φ ( ) (r , r ) i , j <k 3 ijk ij ik (7) ⎛σ + σ j φij(2 ) (rij ) = A⎜⎜ i ⎝ rij ηij ⎞ ⎟ ⎟ ⎠ + Zi Z j rij e − rij λ − α i Z 2j + α j Z i2 4 ij r e − rij ξ r r ⎞ ⎞⎛ rij ⋅ rik ⎛ ξ r r ξ ⎟⎜ φijk (rij , rik ) = Bijk exp⎜⎜ − cos θ ijk ⎟Θ(rij − r0 )Θ(rik − r0 ) + ⎟ ⎟⎜ ⎠ ⎝ rij − ro rik − r0 ⎠⎝ rij rik (3 ) (8) (9) 式(7)の 2 粒子間ポテンシャルは,式(8)の 2 体間相互作用項と式(9)の 3 体間相互作用項からなる.式(9) の Θ(rij-r0), Θ(rik-r0)はカットオフ半径 r0 に関する階段関数である.a-SiO2 のモデル化の計算に用いた各パ ラメータは Vashishta らの文献[5]にある値を用いた.以上より,式(1)を保存量として式(2)から式(4)の運 動方程式を解くことでアモルファス構造のモデル化を行った.初期構造として α-cristobalite 型の SiO2 結 晶を 1000K の温度でアニール後,約 200K の低温まで急冷することにより a-SiO2 の準安定構造をモデル 化することを試みた. また, SiO2 結晶中の酸素欠陥の挙動についても同様のシミュレーションを行った. 古典 MD 計算においては 104 オーダの原子を扱うことが可能であるが,続けて第一原理計算をすること 数時間発展の数値積分にはGear のPredictor-Corrector を考慮してSi32O64 の96 原子の系を扱うにとどめた. -18 を用い,そのときのタイムステップは 2.5×10 s とした.分子動力学計算のフローチャートを図 1 に示 す. 2.2 第一原理計算 次に,分子動力学計算で得られた構 造をさらに緩和するために,第一原理電子状態計算を用 いた.具体的には周期境界条件を課した系に対して密度 汎関数法によるバンド計算を行った.その際,交換相関 汎関数には Perdew と Wang により定式化された一般化 勾配近似(GGA)[6,7]を採用した.電子の波動関数は平 面波で展開されエネルギーに換算して 50Ry までをそ の基底とした.電子密度の計算においては,Γ点にお 図 1 MD 計算のフローチャート ける電子状態計算により運動量空間の状態を近似して いる.本計算で用いた α-cristobalite 型 SiO2 の超格子は 9.956Å×9.956Å×13.896Å の体積を持ち,これは基本単位格子 8 倍の体積であることからΓ点による近似 は妥当であると考えられる.また,計算の効率化のために,原子の内殻を占める電子についてはポテン シャルとして扱う擬ポテンシャル法を用いている.本計算では Hamann によるノルム保存型擬ポテンシ ャル[8]を採用している.以上のような計算により,古典 MD で得られた構造を第一原理バンド計算によ りさらに緩和した.続いて,このようにモデル化された a-SiO2 の電子状態を詳しく調べるために状態密 度(DOS)の計算を行った.DOS の計算は Lehmann と Taut による線形補間の方法[9]に基づき,運動量空 間を連続的に占める k 点の状態を近似した. 結晶とアモルファス構造の DOS の違いを比較することによ り,構造の中の結合状態の変化を明らかにできる.a-SiO2 では,結晶とは異なり,ダングリングボンド を持つ Si および O 原子が存在することが予想され,それらの影響が DOS に現れるはずである. 2.3 量子エネルギー密度 バンド計算の結果から系全体の示す電子状態が明らかになるが,さらに各点 における電子状態を調べるために, 運動エネルギー密度, 張力密度, およびストレステンソル密度[10,11] を定義する.これらの量をみることにより,各点に分布する電子密度の示す物理的・化学的特性を明ら かにすることができる.運動エネルギー密度,張力密度,およびストレステンソル密度はそれぞれ,式 (10),式(11),および式(12)のように表わされる. [ ] r r r r r h2 nT (r ) = − ν i ϕ i* (r )Δϕ i (r ) + Δϕ i* (r )ϕ i (r ) ∑ 4m i (10) r r r r ⎡ * r ∂Δϕ i (r ) ∂ϕ i* (r ) r ∂Δϕ i* (r ) r r h2 * r ∂ϕ i (r ) ⎤ τ (r ) = ∑ν i ⎢ϕi (r ) ∂x k − ∂x k Δϕi (r ) + ∂x k ϕi (r ) − Δϕi (r ) ∂x k ⎥ (11) 4m i ⎣ ⎦ Sk r τ Skl (r ) = r r r r r r ⎡ * r ∂ 2ϕ i (r ) ∂ϕ i* (r ) ∂ϕ i (r ) ∂ 2ϕ i* (r ) r ∂ϕ i* (r ) ∂ϕ i (r ) ⎤ h2 ( ) ( ) ν ϕ ϕ r r − + − ∑ i ⎢ i ∂x k ∂x l ∂x k ∂x l ∂x k ∂x l i ⎥ 4m i ⎣ ∂x l ∂x k ⎦ (12) ここで,{k,l}={x,y,z}である.ストレステンソル密度はエネルギー密度の次元を持つが,単位面積あたり の力とみることもできることから,各点の電子密度に及ぼされる引張りおよび圧縮の応力を意味する. また,式(11)の張力密度は式(12)の発散をとることから求めることができ,各点での電子密度にかかる力 である.式(11),(12)の量は電子の運動に注目しているところから求まる量であるが,平衡状態において はローレンツ力とつりあっている. 3. 結果と考察 3.1 MD 計算による a-SiO2 と酸素欠陥の挙動について前述の 計算方法による MD 計算の結果の一例を示す.α-cristobalite 型 SiO2 の超格子 Si32O64 の組成を持つ構造を 300K から 1000K まで熱した後,急冷した結果得られた準安定の a-SiO2 の構造 を図 2(a)に示す.この構造における結合長および結合角の平 図 2 (a)Si32O64 と(b)Si32O61 の MD 計算 均値はそれぞれ Si–O が 1.618Å,∠OSiO が 109.5°となり と第一原理計算により得られた準安定 Vashishta らの結果とよく一致するものであり,また結晶の値 構造 と大差はない.しかし,∠SiOSi は 169.5°となり,そのため 密度が 1.95g/cm3 となり実測から知られている密度 2.20g/cm3 に比べて小さな値が得られた.この結果は式(7)-式(9)に用い られるポテンシャルのパラメータに依存するところが大きい が,本研究で行った MD 計算の方法が Vashishta らの方法と異 なることにもよる. Vashishta らは a-SiO2 の密度が 2.20g/cm3 [5] 図 3 酸素欠陥のパーコレーションの模 であること,すなわち体積一定であることを束縛条件とした 式図 MD 計算によりアモルファス構造の示す実測データにパラメ ータを最適化しているが,我われの MD 計算では温度と圧力を束縛条件とした.実際の結晶は温度を上 げるに従ってその格子振動が励起され,次第に結晶構造に乱れが生じることから体積膨張を起こすのが 一般的である.その意味では結晶のアモルファス構造を求めるために MD 計算を行う場合は,圧力と温 度の制御のもとで計算を行うべきである.以上より,本計算で用いた経験的ポテンシャルは,我われの 方法で計算する場合にはさらなる最適化が必要であるといえる.図 2(b)に Si32O61 の組成に対して同様の 計算を行った結果を示す.この結果においてもやはり,∠SiOSi が 167.1°と大きくなり,密度が小さく なる傾向にあった.Si32O61 の組成について,初期構造において1個の Si 原子に 2 個の O 原子の欠陥が ある構造から MD 計算を行った結果,図 3 の模式図で示すように隣り合う Si 原子間で O 原子が移動し ているということがわかった.この構造は第一原理計算によって緩和した後にも安定に存在する構造で あることから,熱励起により活性化された O 原子またはその欠陥はパーコレーションを起こすことが示 された. 3.2 a-SiO2 の電子状態 次に,古典 MD 計算と第一原理計算の結果得られた SiO2 の準安定構造について の DOS の計算結果を図 4 に示す.図 4 には Si32O64,Si32O63,Si32O61,および Si32O59 の組成について計 算した結果を示す.Si32O64 の組成は SiO2 の組成となっているが,結合長と結合角についての乱れのため に最安定な結晶に比べてピークに幅が見られる.また,酸素欠陥を含む場合の DOS においては価電子 帯と伝導帯に SiO2 では見られなかった準位が現れた.これは O 原子の欠陥のために Si 原子上に残され たダングリングボンドに起因するものである.しかし,その準位の位置はさまざまであり,Si32O64(SiO2) の組成から O 原子を 1 個抜き取った Si32O63 の構造では価電子帯上端に準位が現れ,Si32O61 の構造では 伝導帯下端に,さらに Si32O59 の場合はバンドギャップ間に準位が現れる.これらの違いは O 原子を抜き 取ることによりモデル化した酸素欠陥を MD 計算により緩和した結果であり,組成による欠陥の安定構 図4 (a)Si32O64,(b)Si32O63,(c)Si32O61,および(d)Si32O59 の構造の示す DOS.(a)の VT,CB はそれぞれ価電 子帯上端および伝導帯下端を表す. 造に違いが現れたといえる.Si32O63 の構造のように価電子帯上端に準位が現れるのは酸素欠陥により生 じたダングリングボンドが周囲の原子により消されていると考えられる.一方で,ダングリングボンド が緩和されずに残る場合は,Si32O59 の DOS に見られるようにバンドギャップ間に不純物準位のように 出る.ダングリングボンドが消される場合でも反結合性の不安定な結合の生成による場合は,その準位 は伝導帯側に現れると考えられ Si32O61 の DOS で示されるようになる.以上のように,本研究では圧力 と温度を考慮した MD 計算と第一原理計算により,a-SiO2 の中の電子構造をモデル中に再現することを 試み,その結果,実際に存在すると考えられている構造の電子状態を見ることを明らかにした. 3.3 量子エネルギー密度 上述の MD 計算と第一原理計算により得られた電子状態についての DOS 計算 の結果をさらに局所的に考察するために式(10),式(11),および式(12)に示されているように,運動エネ ルギー密度,張力密度,およびストレステンソル密度の計算を行った.図 5 に Si32O63,Si32O61 の組成に 見られた特徴的な構造を示す.図 5(a),(b)は,Si32O63 に見られた価電子帯上端の準位に対応する Si–Si 間の結合を示す運動エネルギー密度とストレステンソル密度である.これは O 原子を抜き取った後, MD 計算により緩和したところ生じたものである.図 5(a)の運動エネルギー密度を見ると,Si 原子間に 密度があることから 2 原子間の電子の交換を読み取ることができ,2 原子の境界が張力密度で示されて いる.その境界面に密度が存在することで,一方の原子に束縛されている電子が他方に乗り移ることが 明らかにされている. 図 5(b)のストレステンソル密度を見ると Si–Si 間に正の値の大きな領域が存在し, 各点での固有ベクトルは Si–Si 間で連続的であると読み取れることから Si 原子間の共有結合が示唆され る.同様の図を Si32O61 の系について図示したものが図 5(c),(d)である.この場合の DOS は伝導帯下端 に準位が現れたもので,その準位に対応するものを図示したものである.これらの図では Si 原子間の距 離が 4.5Å と長いために強い結合は考えられないが,Si-Si 間の共有結合的は結合を見ることができる. この結合距離と結合の弱さが図 4 に示された DOS に反映されており,結合の不安定さのために価電子 帯下端に準位が生じたものと考えられる. 図 5 Si32O63 の Si–Si 間に見られる(a)運動エネルギー密度と張力密度および(b)ストレステンソル密度. Si32O61 の Si–Si 間に見られる(c)運動エネルギー密度と張力密度および(d)ストレステンソル密度. 4. まとめ 本研究においては,圧力と温度を考慮した分子動力学計算と第一原理計算による電子構造の解析により, アモルファス構造のモデル化についての方法を述べた.その背景にはゲート絶縁膜として使われている a-SiO2 の構造と電子状態の解明という課題がある.今回は,a-SiO2 のモデル化と絶縁破壊の原因のひとつに 挙げられる酸素欠陥の挙動に注目して計算を行った.その結果,アモルファスとしての準安定構造の特徴と酸 素欠陥の影響について議論することができ,古典的 MD 計算と第一原理計算を組み合わせることで a-SiO2 の 絶縁膜のモデル化に対して向上が得られたと考える.また,状態密度という系全体の電子構造から,エネルギ ー密度といった局所的な構造についての焦点を絞った議論を行うことにより,実測で得られる平均化したデー タに対しての局所的な原因解明について示唆を与えることができることを示した.今後は,題目にも挙げている ように,今回の研究で得られたモデルをもとにして,さらに電場応答の計算を行うことで絶縁破壊の原因を解明 したいと考えている. 参考文献 [1] S. Yoshida, K. Doi, K. Nakamura, and A. Tachibana, Appl. Sur. Sci. 216 (2003) 141. [2] K. Doi, K. Nakamura, and A. Tachibana, Appl. Sur. Sci 216 (2003) 463. [3] S. Nosé, J. Chem. Phys. 81 (1984) 511. [4] M. Parrinello and A. Rahman, J. Appl. Phys. 52 (1981) 7182. [5] A. Nakano, R. K. Kalia, and P. Vashishta, J. Non-Cryst. Solids 171 (1994) 157. [6] J. P. Perdew and Y. Wang, Phys. Rev. B 45 (1992) 13244. [7] J. P. Perdew, J. A. Chevary, S. H. Vosko, K. A. Jackson, M. R. Pederson, D. J. Singh, and C. Fiolhais, Phys. Rev. B 46 (1992) 6671. [8] D. R. Hamann, Phys. Rev. B 40 (1989) 2980. [9] G. Lehmann and M. Taut, Phys. Stat. Sol. (b) 54 (1972) 469. [10] A. Tachibana, J. Chem. Phys. 115 (2001) 3497. [11] A. Tachibana, J. Mol. Model. 11 (2005) 301. 2 次元非線形格子モデルにおけるエネルギー局在構造の解析 工学研究科航空宇宙工学専攻 土井祐介* Abstract: Structure of intrinsic localized modes (ILMs) in two dimensional anharmonic lattice systems is investigated. Numerical solutions of static ILMs are calculated by means of finding the periodic orbits in the phase space using the Newton-Raphson method. We find two types of the localized structure: (a) ILM with quasi-one dimensional symmetry and (b) ILM with two dimensional symmetry. Shape of the ILMs is affected by the microscopic structure and shape of interaction potential. Key words: Intrinsic localized mode, Two dimensional lattice, Energy Localization 1. はじめに 非線形格子モデルにおいて出現する非線形局在モード (Intrinsic localized mode/ILM)[1] は,非可 積分系において出現することから理学・工学における様々な系での出現が期待されている.実際 に,光導波路アレイでの光の局在現象 [2],マイクロカンチレバーアレイの局在振動現象 [3] など の実験において ILM が出現することが報告されている. ILM と材料の関係について考えると,結晶構造は原子・分子の周期構造であり,格子モデルと してモデル化することが可能である.またそれらの相互作用は本質的に非線形性を有しており,非 線形性が強く顕在化する状況においては材料中に ILM が出現することが期待される.したがって, 実際の材料中のダイナミクスにおいて ILM の果たす役割を調べることは重要である. 一方で,従来 ILM の研究の多くは 1 次元格子モデルにおいて解析がなされてきている.しかし ながら,ILM の生成および維持のメカニズムは系の次元数には依存せず,2 次元,3 次元系におい ても ILM が存在しうることが示されている.実際に 2 次元非線形格子モデルにおいていくつかの ILM の生成・維持のシミュレーションが行われている.Marı́n らは相互作用として相互作用ポテ ンシャルが Lennard-Jones ポテンシャルで与えられる 3 角形格子系において,数個の格子点の擾乱 から長時間安定に存在する移動型の ILM が形成されることを示した [4].また,池田は相互作用 が Morse ポテンシャルである 3 角形および 4 角形格子系において最も高い振動数をもつフォノン モード (Zone Boundary Mode/ ZBM) の変調不安定から 2 つのタイプの局在構造が出現することを 示している [5].さらに,山寄らは分子動力学シミュレーションによってグラフェンシートにおい て ILM が出現しうることを示した [6]. 以上のように 2 次元格子系において比較的長時間にわたって安定に存在する局在が出現するこ とが数値的に確かめられている.その一方で,理想的な格子モデルにおいても高精度な ILM の数 値解は求められていない.このような解を求めることによって,ILM の精密な構造の解析が可能 になり,上記のシミュレーション結果の理解に寄与すると考えられる.そこで本研究では 2 次元 格子モデルにおける粒子の面内振動に着目して ILM の高精度な数値解を計算し,構造を解析する. 2. モデル 本研究では,図 1 に示す 2 次元正方格子モデルを考える.格子点 は x 方向,y 方向にそれぞれ等間隔に周期 d = 1 で配列している. x 方向に i 番目,y 方向に j 番目の格子点を (i, j) 格子点とする.そ れぞれの格子点の釣り合いの位置からの変位は ui,j = (xi,j , yi,j ) とする.各格子点は最近接格子点 4 点および第 2 近接格子点 4 点と相互作用をすると考える.格子間の相互作用ポテンシャル としては Fermi-Pasta-Ulam (FPU) β 型とする.この場合,系の _______ * 平成 17 年 7 月より大阪大学大学院 工学研究科 知能・機能創成工学専攻 助手 (i-1, j+1) (i, j+1) (i+1, j+1) y (i-1, j) (i, j) (i+1, j) (i-1, j-1) (i, j-1) (i+1, j-1) x 図 1: 2 次元格子モデル Hamiltonian は H = X 1X 2 u̇i,j + [V (|ui,j − ui+1,j | − L1 ) + V (|ui,j − ui,j+1 | − L1 ) 2 i,j i,j + αV (|ui,j − ui+1,j+1 | − L2 ) + αV (|ui,j − ui−1,j+1 | − L2 )] , 1 2 1 4 V (r) = r + βr , 2 4 (1) (2) で与えられる.ここで L1 , L2 は最近接格子,第 2 近接格子の釣り合いの長さ,α は第 2 近接相互 作用の大きさを表すパラメータである.以下では L1 = 1 とする. 式 (1) から各格子点の運動方程式は d2 ui,j dt2 = F(ui,j , ui+1,j , 1) + F(ui,j , ui−1,j , 1) + F(ui,j , ui,j+1 , 1) + F(ui,j , ui,j−1 , 1) + αF(ui,j , ui+1,j+1 , L2 ) + αF(ui,j , ui−1,j+1 , L2 ) + αF(ui,j , ui+1,j−1 , L2 ) + αF(ui,j , ui−1,j−1 , L2 ), (3) となる.F(u1 , u2 , d) は格子間に作用する力である. u2 − u1 , |u2 − u1 | f (r) = r + βr3 . F(u1 , u2 , d) = −f (|u2 − u1 | − d) (4) (5) この系の線形分散関係は x, y 方向の波数をそれぞれ k, l として ω 2 = −4A cos k cos l + C − (cos k + cos l) ± q (cos k − cos l)2 + 16B 2 sin2 k sin2 l, (6) で与えられる.ただし, √ √ √ α (4 − 2L2 + 16β − 18 2L2 β + 12L22 β − 2L32 β), 4 √ √ α √ ( 2L2 + 8β − 6 2L2 β + 2L32 β), B = 4 √ √ √ C = (2 + 4α − 2L2 α + 16αβ − 18 2L2 αβ + 12L22 αβ − 2L32 αβ), A = (7) (8) (9) である.この式によって得られる分散曲線の例を図 2(a) に示す.式 (6) は k = π, l = 0 または k = 0, l = π の時に最大値 2 ωmax = C + 2 − 4A, (10) をとる.このときの系の変位パターンを図 2(b) に示す. 4.5 4 3.5 3 ω 2.5 2 B 1.5 1 O 0.5 0 O A B (a) 線形分散関係 O A k = π, l = 0 k = 0, l = π (b)ZBM の変位パターン 図 2: 2 次元正方格子モデルの線形分散関係および ZBM の変位パターン 3. 数値解析 ILM は空間的に局在した周期解であり,位相空間上では周期軌道をとる.本研究では式 (1) で記 述される系において線形定在波の最大角振動数 ωmax よりも大きな角振動数 ωILM をもつ周期軌道 を探す.式 (1) で規定される系において時刻 0 に位相空間内の点 X0 の時間発展は X(t) = F(t, X0 ), (11) とあらわす.X0 が内部振動数 ωILM の ILM のある瞬間の位相空間上での位置とする,この点は周 期軌道上にある.したがって n を整数として 2π F(t, X0 ) = F(t + n , X0 ), (12) ωILM である.いま,位相空間上の点 X が式 (6) に従って,時刻 2π ωILM 後に移動する点 X0 への写像を X0 = MωILM (X), (13) MωILM (X) − X = 0, (14) とするならば,式 (12) の関係は となる.これは X に対する方程式であり,これを満たす X を求めることが内部振動数 ωILM の ILM を求めることに対応する.ここでは式 (14) の解を Newton-Raphson 法を用いて計算する.反 復計算中には写像 MωILM の計算を行う必要があるが,MωILM は解析的な形では与えられず,式 (3) を積分することによって求める.数値積分には 6 次のシンプレクティック数値積分を用いる. Newton-Raphson 法の初期解としては,ILM の局在性が非常に強い場合(3 格子点または 4 格子 点に局在が集中している場合)の近似方程式を求め,その解を用いている. 4. 結果および考察 3 3 3 2 2 2 1 1 1 0 -1 -1 -1 -2 -2 -2 -3 -3 -3 -2 -1 0 x 1 2 3 4 1.5 j=1 1 0.5 0 -0.5 -20 -3 -2 -1 0 x 1 2 3 4 -4 1.5 Displacement xi,j 1 -3 -4 j=0 -15 -10 -5 0 i 5 (a)α = 0 10 15 20 j=1 1 0.5 0 -0.5 -20 -3 -2 -1 0 x 1 2 3 4 5 10 15 20 1.5 Displacement xi,j -4 Displacement xi,j 0 y 0 y y 4.1 準 1 次元 ILM 図 3 に準 1 次元的な ILM の計算結果を示す.これらの図から格子点の x 軸方 向の振動が y 方向に比べて卓越していることが分かる.また x 方向の振動も一つの列 (j = 0) に集 中していることが分かる.図 3(a) は第 2 近接相互作用がない場合の解である.この場合は 1 次元 系と同様に左右対称な局在構造が出現する.一方,第 2 近接相互作用を考慮した場合,図 3(b)(c) のように局在構造の対称性が失われることが分かる. j=0 -15 0.5 0 -10 -5 0 i 5 10 15 (b)α = 1, L2 = 1, ω = 2.4 図 3: 準 1 次元 ILM の構造 20 j=1 -0.5 -20 j=0 -15 -10 -5 0 i (c)α = 1, L2 = 1, ω = 4.2 0.8 d2=√2−0.1 d2=√2−0.2 d2=√2−0.3 d2=1.0 0.7 0.6 0.5 0.4 0.3 0.2 0.1 3 3.5 4 Angular frequency ωILM 4.5 図 4: 準 1 次元 ILM の内部振動数 4.2 2 次元 ILM 数値計算によって,2 次元的な対称性をもつ ILM も存在することを確かめた.その例を図 5 に示す.これら の ILM は x, y 方向に同じオーダの振幅で振動する.非常に局在 性が強い場合は正方格子を構成する 4 点が正方格子の中心に向 かって振動する形をとる.この 2 次元 ILM についても,線形の 固有モードの振動パターンの観点から関連する固有モードが存 在すると考えられる. 3 2 1 y Amplitude for x-direction 図 4 に系に励起される準 1 次元 ILM の振幅と内部振 動数の関係を示す.系のパラメータ L2 の大きさが異な る 4 つの場合について示す.図中の直線はそれぞれの条 件における線形分散固有角振動数の最大値である.それ ぞれの結果から励起された ILM は線形の最大固有角振動 数モード (ZBM) よりも大きな内部振動を持つことが分か る.ZBM の振動パターンは図 2(b) に示している形であ り,準 1 次元 ILM の局在構造はこの振動パターンに由来 するものと予想される.また,L2 が小さくなるにつれて 生成される ILM の内部振動数が大きくなる傾向があるこ とが分かる.第 2 近接相互作用は系に 2 次元性を導入す ることによって大きく現れる効果であることから,2 次 元性の影響によって ILM の振動数が増大するといえる. 0 -1 -2 5. おわりに -3 -2 -1 0 1 2 3 4 1 2 3 4 x 3 2 1 y 本研究では,2 次元非線形格子系における ILM の構造を数値 的に解析した.その結果,線形の ZBM に関連する準 1 次元的 な ILM が存在することが明らかになった.また,準 1 次元的な ILM は 1 次元格子における ILM とよく似た構造,励起メカニズ ムであるものの,その局在構造および内部振動数は系の 2 次元 性の影響をうけて大きく変化することを確認した.さらに,2 次 元的な対称性を持つ 2 次元 ILM が出現しうることを確かめた. これらの構造については,解として存在しうるが,その安定 性,および励起のメカニズムについての詳細は明らかではない. 今後はこれらの問題を解明していく予定である. 0 -1 -2 -3 -2 -1 0 x 図 5: 2 次元 ILM の構造 参考文献 [1] A.J Sievers and S. Takeno, Intrinsic localized modes in anharmonic crystals, Phys. Rev. Lett., 61 (1988), 970. [2] H.S. Eisenberg et al, Discrete spatial optical solitons in waveguide arrays, Phys. Rev. Lett, Vol.81 (1998), 3383. [3] M. Sato et al, Observation of locked intrinsic localized vibrational modes in a micromechanical oscillator array, Phys. Rev. Lett., Vol.90 (2003), 044102. [4] J.L. Marı́n et al, Localized moving breathers in a 2D hexagonal lattice, Phys. Lett. A, 248 (1998), 225. [5] 池田, Morse potential を持つ格子系における非線形局在モードの生成, 京都大学数理解析研究 所講究録 1368 (2004), 189. [6] 山寄他, 分子動力学法によるグラフェンシートにおける非線形局在モードの解析, 第 18 回計 算力学講演会講演講演論文集 (2005), 657. 材料における原子領域のトポロジー 工学研究科マイクロエンジニアリング専攻 中村 康一 Abstract: The partitioning a material system into atomic regions in terms of electronic-state theory has been examined. In addition to the conventional atoms-in-molecules (AIM) method by Bader based on electron density, new methods of regional partitioning were introduced based on quantum energy density such as the electronic kinetic energy density and tension density. As a preliminary work to application to materials, we have applied the regional partitioning to HeH+ molecule and HeLi+ diatomic system by ab initio electronic-state calculation. We have discussed topological dependency of atomic regions on basis sets and exchange-correlation functionals for each method of partitioning, and computed amount of charge transfer between regions, that is related to atomic charge of each region. Furthermore, application to LiF crystal system in lattice vibration has been carried out by means of first-principle calculation with periodic-boundary LCAO wave functions. Topology of atomic regions will be characterized by physical meanings due to the base for partitioning in the system. Key words: Regional partitioning, First-principle calculation, Atomic charge, Basis set, Exchange-correlation functional, Vibration 1. はじめに 材料を原子の集合体として見た場合、その材料の物性は数多くの「原子自らがもつ性質」と「原子間 相互作用」から成り立っていると解釈できる。系が原子の集合体であるとき、もはや原子は「点」では ありえず、有限または無限の体積をもった領域を占めることになり、原子自らがもつ性質と原子間相互 作用の本質がそれぞれ領域内の相互作用、領域間の相互作用(とりわけ領域境界における性質)として 与えられる。原子が占める領域を決める方法として、電子密度を基準とするトポロジカルな電子状態理 論手法が Bader により提案されており(Atoms-in-Molecules (AIM)法[1]) 、分子系における一般的な方法 として広く用いられているが、原子電荷や励起状態でのトポロジーに定性的・定量的な不具合がみられ るほか、材料系への応用はあまり進んでいないのが現状である。本研究では AIM 法のほか、我々のグル ープで開発したプログラム[2,3]とリンクさせて量子エネルギー密度[4]を基準とした領域分割を行い、波 動関数展開方法等の計算条件依存も含めて材料における原子領域のトポロジーと含有する電子物性を解 析することを目的とする。 2. 領域分割の手法 2.1 AIM理論 古くから、全系のエネルギーを原子や部位ごとに割り当てる、いわゆる「エネルギー分 割」という概念は、化学的直感に対応する重要な解析手法として多くの研究者により議論され、今日で もさまざまな方法が提案されつづけている。AIM理論[1]は、全系を領域分割する手法としてシンプルか つ整然であるとともに、原子領域や化学結合など分子における多くの基礎的概念を踏襲する手法として 広く認知されており、Gaussianシリーズ[5]などの著名なプログラムパッケージに組み込まれているほか、 より詳細な情報を得るためのパッケージもいくつか発表されている[6,7]。AIM法の最大のポイントは観 測量である電子密度のトポロジカルな情報を用いることであり、量子化学計算によって得られる電子密 度を用いることによって領域分割の境界や分割に基づく領域物性量がユニークに決定される。 量子化学計算による電子密度ρ(r)は、一般に自然軌道ψi(r)とその占有数νiを用いて ρ (r ) = ∑ν iψ i* (r )ψ i (r ) (1) i のように定義できる。通常、ρ(r)は原子核のある位置のみで極大となり、原子核から離れるにつれ急速 に減少していく。電子密度の勾配∇ρ(r)は電子密度の増加する方向を向き、空間上のある位置から∇ρ(r) に沿って空間を動いていくと(密度勾配線) 、その密度勾配線は原子核上にあるいずれかのρ(r)の極大点 (アトラクター)に到着することになる。原子核A上のアトラクターに到着する密度勾配線の始点の集 合が原子Aに属する原子領域として定義される。 原子領域どうしの境界面S上の点rSにおいては∇ρ(rS) · n(rS) = 0が成り立ち(ただし、n(rS)は境界面Sの 法線ベクトル) 、また∇ρ(r) = 0となる位置を臨界点(critical point)と呼ぶ。領域分割の手順としては、隣 接する2つの原子領域のアトラクターを結ぶ密度勾配線上の臨界点(bond critical point)を探し、その位 置から–∇ρ(r)の向きに沿ってρ(r)の極小点をめざして空間を動く方法が基本となる。 AIM法は身近な電子密度に基づくことできわめて直感的にわかりやすい手法であるが、この方法に従 って得られた原子領域内の電子密度を積分するという形で導出した原子電荷(atomic charge)が実際に想 定されるatomic chargeと大きくかけ離れる場合があるなど、トポロジーに定性的・定量的な不具合がみら れることがある。領域分割においては電子密度に基づかなければならない必然性はなく、別の物性量に 基づいて行う領域分割も議論する余地がある。 2.2 量子エネルギー密度 領域密度汎関数理論に基づき、電子の運動や電子が関与する相互作用を記述 する局所物性量として、運動エネルギー密度 nT(r)が式(2)のように自然軌道を用いて定義できる[4]。 nT (r ) = ⎛ ⎧ h2 ⎫ ⎧ h2 ⎫⎞ 1 * * ⎜ ν ψ r ψ r ψ r ( ) ( ) ( ) − Δ + Δψ i (r ) ⎬ ⎟⎟ ⎨ ⎬ ⎨− ∑ i⎜ i i i 2 i ⎝ ⎩ 2m ⎭ ⎩ 2m ⎭⎠ (2) 系空間は nT(r)の正負によって、電子が古典力学的に自由に運動可能である領域 RD (electronic drop region; nT(r) > 0)と量子力学的なトンネル効果を用いてのみ運動を許される領域 RA (electronic atmosphere region; nT(r) < 0)に分割でき、2原子間を連結する RD の存在によって化学結合が表現されるなど、電子間相互作 用の概要が視覚化される。この物性量に基づいて領域分割を行うと、導出される領域物性量に電子運動 や相互作用の影響が暗に含まれることが期待できる。領域分割は AIM 法と同様にして、原子核 A 上の アトラクター(この場合、nT(r)の極大点)に到着する運動エネルギー密度勾配線の始点の集合を原子 A に属する原子領域として定義する。 また、系にはたらく電子的な総ての力の密度 FS(r)は電子テンション密度τS(r)と電子外力密度 XS(r)の 和であり、とりわけτS(r)は純粋に量子力学の由来からその成分が τ Sk (r ) = ⎛ * ∂Δψ i (r ) ∂ψ i* (r ) ∂Δψ i* (r ) ∂ψ (r ) ⎞ h2 ⎜ − Δ + ( ) ( ) ν ψ r ψ r ψ i (r ) − Δψ i* (r ) i k ⎟⎟ ∑ i⎜ i i k k k ∂x ∂x ∂x ∂x ⎠ 4m i ⎝ (3) と導かれ、定常状態において FS(r) = 0 となることから電場 E(r)によって電子にはたらく力と釣り合って いる[4]。テンション密度に基づく電子にはたらく力をベースにした領域分割も興味が持たれるところで あり、この場合、–τS(r)の向きに沿ってに沿って空間を動いた軌跡が原子核 A 上に到着するときの軌跡 始点の集合を原子 A に属する原子領域として定義する。 3. 応用計算 3.1 HeH+分子 材料系について領域分割の応用計算を行う前に、はじめに2原子分子系についての計算 を行い、基本事項を確認した[8]。HeH+分子の Hartree-Fock (HF)電子状態についての領域分割結果を図 1 に示す。HeH+分子は最も簡単な異核2電子系であり、化学結合による安定構造が存在するとともに、原 子間距離が相当長い場合でも HF 法や密度汎関数法で用いられる単一行列式による波動関数に信頼性が ある等、きわめて解析に適した系である。 HeH+分子の安定構造において、AIM 法による領域分割では HF/6-311++G レベルで辛うじて He–H 原 子核間に bond critical point が現れて He 領域と H 領域に分割できるものの、HF/STO-3G レベルや HF/6-311++G(3df,3pd)レベルでは bond critical point が現れない。一方、nT(r)を基準とする領域分割法(以 下、nT(r)法と略す)においては、He 領域と H 領域がともに閉領域として現れ、それぞれの原子核から 離れた地点にどの原子領域にも属さない空間が存在する。また、τS(r)を基準とする領域分割法(以下、 τS(r)法と略す)においてはそのトポロジーにきわめて大きな基底関数依存性が見られた。 3 3 (B) (C) 2 3 (C) 1 1 (C) (C) (B) He 0 -1 (C) H (C) -1 (B) -2 2 (B) He 0 (C) (A) 2 1 (C) (B) H -1 (B) -2 -1 0 1 2 3 4 5 -3 -2 -1 R(He–H) = 0.929 Å 0 1 2 3 4 -3 5 (C) 1 (B) (A) He 0 -1 (A) 2 (C) (B) 1 H He (B) H (B) -2 -2 -2 -3 -3 -1 0 1 2 3 R(He–H) = 1.500 Å 4 5 6 He 0 -1 -2 1 2 3 4 5 (A) 1 0 (B) 0 (C) 2 -1 -3 -1 3 (C) (C) -2 R(He–H) = 0.771 Å 3 2 (B) R(He–H) = 0.773 Å 3 (C) -3 -3 -3 (B) H -2 -2 -3 (B) He 0 (C) (B) H (B) -3 -3 -2 -1 0 1 2 3 R(He–H) = 1.500 Å 4 5 6 -3 -2 -1 0 1 2 3 4 5 6 R(He–H) = 1.500 Å (a) (b) (c) 図 1 最安定構造および He–H 距離が 1.5 Å における HeH+分子の原子領域分割:(a) HF/STO-3G レベル、 (b) HF/6-311++G レベル、(c) HF/6-311++G(3df,3pd)レベル。以降の図も含めて、境界面(A)は AIM 法、(B) は nT(r)法、(C)はτS(r)法によるものであり、背景の黒またはグレーの領域は RD、矢印はτS(r)の向きをそれぞ れ表す。座標軸の数値単位は bohr。 HeH+分子系では電子の大半が He 原子核付近に存在するため、ρ(r)の基底関数依存性は H 原子核付近 でほとんど見られないが、nT(r)やτS(r)は図 1 に示したように H 原子核付近で大きな基底関数依存性が見 られる。これは nT(r)やτS(r)が自然軌道の2階微分を含むため、基底関数が微分されることでその指数パ ラメータが乗じられるという影響が露になっている。nT(r)やτS(r)の変化に応じて、これらそれぞれを基 にした原子領域の境界も大きく変化する。 実際の一般的な物性解析においては、 外殻軌道のスプリット、 さらには分極関数や diffuse 関数の追加なしには定量的な議論ができず、原子領域について大きな基底関 数依存性が見られることは自然な現象といえる。 3.2 HeLi+系 通常、材料系の第一原理計算は密度汎関数法に基づく計算となる。密度汎関数法は比較 的少ない時間とメモリのコストで実測とよく対応する計算結果を導くが、最大の問題点の1つが長距離 での電子間相互作用が満足に記述できないことであり、長距離相互作用の精度については HF 法よりも 劣るためしばしば HF 法の交換エネルギーを混成した交換相関汎関数が用いられる[9]。原子間距離が相 当長い場合でも単一行列式による波動関数に信頼性がある条件で長距離電子間相互作用が議論できる最 も簡単な系が HeLi+系である。HF 法および密度汎関数法を用いた HeLi+安定構造についての領域分割の 例を図 2 に示す。また、本研究で検討した密度汎関数法のリストと 6-311++G(d,p)基底関数を用いた安定 構造での He–Li 距離、各領域分割方法における He–Li 原子核間上の critical point と He 原子核との距離を 表 1 に示す。 HeLi+安定構造での He–Li 距離は、交換項を HF 法の交換エネルギーから Slater 交換汎関数に変えるこ とによりわずかに短くなり、さらに Becke の勾配補正を加えると逆にはっきりと長くなる結果が得られ た。また、相関汎関数(とりわけ LYP 相関汎関数)は He–Li 距離を短くなる方向に寄与することも示さ れた。このように用いる汎関数によって He–Li 距離には 0.2 Å 以上の違いが生じ、密度汎関数法による 長距離相互作用記述の難しさを表しているが、いずれにせよ He–Li 距離はおよそ 2 Å もあり、図 2 に示 したように He および Li 原子核上にあるそれぞれの RD は互いに連結しておらず、nT(r)の立場では化学 結合をしているとはいい難い。これに対応して、nT(r)法においては、比較的 He–Li 距離の長い HFB 法、 BVWN 法、BP86 法、および BPW91 法の計算結果に関してそれぞれが閉領域である He 領域と Li 領域 3 3 (C) (A) (C) 2 3 (C) (A) (C) 1 (B) 0 He (B) 1 Li 0 2 (B) He (B) (B) 1 Li 0 -1 -1 -1 -2 -2 -2 -3 -3 -3 -2 -1 0 1 2 3 4 5 6 7 -2 -1 0 1 2 3 4 5 6 7 -3 (B) 1 He (B) 1 Li (B) 0 He Li (B) -2 -2 -3 -3 1 2 (d) 1 3 4 5 6 7 2 3 4 5 6 7 5 6 7 (C) (A) (C) (B) (B) 0 -2 0 0 1 -1 -1 -1 2 -1 -2 Li 3 -1 -3 -2 (C) (A) (C) 2 (B) He (c) 3 (C) 2 0 (B) (b) (A) (C) (B) -3 -3 (a) 3 (C) (A) (C) 2 He Li (B) -3 -3 -2 -1 0 1 2 (e) 3 4 5 6 7 -3 -2 -1 0 1 2 3 4 (f) 図 2 密度汎関数法計算の最安定構造における HeLi+系の原子領域分割:(a) Hartree-Fock 法、(b) HFS 法、 (c) HFB 法、(d) BVWN 法、(e) BLYP 法、(f) B3LYP 法。基底関数は 6-311++G(d,p)。 表 1 6-311++G(d,p)基底を用いた密度汎関数法計算による HeLi+最安定構造の He–Li 距離および各領域 分割法での He–Li 原子核間上の He 原子核–critical point 間距離 Optimized Distance between He and critical point Functional Exchange Correlation He–Li length AIM nT(r) τS(r) HF HF none 1.9498 1.0899 1.4417 0.9879 HFS Slater none 1.9028 1.0629 1.3894 0.9559 HFB none 2.0698 1.1662 1.5564, 0.9251 0.9967 Slater+ΔBecke SVWN Slater VWN 1.8536 1.0325 1.3413 0.9465 BVWN VWN 2.0155 1.1314 1.5099, 0.9229 0.9950 Slater+ΔBecke BP86 P86 2.0056 1.1260 1.4985, 0.9325 0.9914 Slater+ΔBecke BPW91 PW91 2.0162 1.1321 1.5093, 1.9249 0.9942 Slater+ΔBecke BLYP LYP 1.9123 1.0693 1.4058 0.9650 Slater+ΔBecke B3LYP 1.8961 1.0590 1.3889 0.9641 Slater, HF, ΔBecke LYP, VWN が互いに隣接しておらず、AIM 法による領域分割では表現できないトポロジーとなる。また、 6-311++G(d,p)基底関数を用いた HeLi+系ではτS(r)の挙動がかなり複雑であるが、τS(r)法による He 領域は シンプルな閉領域で与えられる。一方、He 原子核とは向かい側の Li 原子核から約 1.2 Å の距離の球面 状にτS(r)法による Li 領域とどの原子にも属さない領域との境界がはっきりと観察されたが、調べた系空 間の広さでは境界を確定することができなかった。おそらく He 領域をすべて含むような形で Li 領域が 閉空間として存在すると思われるが、その場合でも He 領域に含まれる分を除いたほぼすべての電子が Li 領域内に存在するので、ここでは単に He 領域の補空間を Li 領域とした。 それぞれの密度汎関数法を用いた HeLi+安定構造での Mullilen atomic charge と各原子領域内の電子密 度を積分して導出した atomic charge を表 2 に示す。各原子領域内の電子密度積分には任意有界領域の領 域電子エネルギー計算プログラム GREG [10]を改良して使用した。Mullilen atomic charge によると He か ら Li+へ 0.051~0.081 の電子移動があるが、AIM 法による領域間にはほとんど電子移動が見られない。 nT(r)法では Mullilen atomic charge によく対応する He 領域からの電子流出が見られたが、この領域分割の 場合はすべての空間をカバーしているわけではなく、さらに議論が必要である。 表2 6-311++G(d,p)基底を用いた密度汎関数法計算による HeLi+最安定構造の Mulliken atomic charge およ び各領域分割法における領域内電子密度積分による atomic charge Atomic charge derived from regional partitioning Functional Mulliken (He/Li) AIM (He/Li) nT(r) (He/Li) τS(r) (He/Li) HF 0.051 / 0.949 0.006 / 0.994 0.060 / 1.058 0.028 / 0.972 HFS 0.079 / 0.921 0.011 / 0.989 0.063 / 1.056 0.033 / 0.967 HFB 0.062 / 0.938 0.008 / 0.992 0.075 / 1.056 0.029 / 0.971 SVWN 0.081 / 0.919 0.011 / 0.989 0.063 / 1.057 0.034 / 0.966 BVWN 0.053 / 0.947 0.006 / 0.994 0.072 / 1.060 0.029 / 0.971 BP86 0.058 / 0.942 0.007 / 0.993 0.073 / 1.059 0.029 / 0.971 BPW91 0.055 / 0.945 0.007 / 0.993 0.073 / 1.059 0.029 / 0.971 BLYP 0.068 / 0.932 0.010 / 0.990 0.061 / 1.059 0.032 / 0.968 B3LYP 0.066 / 0.934 0.009 / 0.991 0.061 / 1.059 0.032 / 0.968 3.3 LiF 結晶の格子振動 材料系への領域分割の適応については、これまで COE 研究で行ってきた LiF 結晶の格子振動を手始めとした[11,12]。格子振動モードは立方晶格子最小セルのモード[11]ではなく、一 般的な単位格子のモード[12]を用いた。周期系の場合、波動関数の展開方法は平面波展開と原子軌道展 開に大きく分けられるが、前者においては原子位置を中心とする同心円状の RD フリンジ構造が見られ るなど[13]フーリエ変換のグリッド数に応じて nT(r)の増減やτS(r)の方向反転が小刻みに生じるため、 nT(r)やτS(r)を基準とする領域分割は実際上不可能となる。したがって、原子軌道展開による自然軌道を 用いて領域分割を検討した。Li 原子に 6-11G 基底、F 原子に 7-311G 基底をそれぞれ使用した場合の領域 分割の結果を図 3 に示す。 Li–F 原子核間において、 AIM 法による bond critical point は Li 原子に近い位置に存在するのに対し、 τS(r) 法においては bond critical point が Li–F 間のほぼ中点に存在する。その結果、AIM 法による領域分割では 凸閉空間の Li 領域を取り囲む F 領域どうしが隣接するが、τS(r)法による領域分割においてはほぼ球状の F 原子領域を取り囲む Li 原子領域どうしが隣接するという完全に逆転したトポロジーの様子が示された。 この系ではほとんど全ての価電子が F 原子付近に存在するため、AIM 法とτS(r)法での領域積分による atomic charge の違いは積分範囲が大きく違うのにも関わらず 0.02 程度しか見られない。価電子の立場か ら系空間を見たとき、もはや Li 原子(正確には Li+イオン)はただのクーロン場でしかなく、τS(r)法に よる領域分割は価電子が詰まったいくつかの孤立した球状の領域をクーロン場の領域が取り囲むという 物理的描像を表現しているとも解釈できる。なお、nT(r)法による原子領域は、Li 原子・F 原子とも大き さが RD よりもわずかに広い球状閉空間であり、bond critical point は見られなかった。すなわち、どの原 子にも所属しない領域が連結して広く分布し、領域間相互作用の議論は困難である。 格子振動の効果について F 領域に着目すると、AIM 法による領域は大きく形が歪むのに対し、τS(r) 法による領域はほとんど変化しなかった。振動に応じてそれぞれの分割法による bond critical point が互 いに接近する等の様子が観察されたが、各領域分割法においてトポロジカルな変化が見られるわけでは なく、基底関数依存等に比べると振動の効果は小さいという結果が得られた。 4. まとめと今後の課題 分子系および周期系モデルについて AIM 法や量子エネルギー密度に基づく原子領域分割を行ったと ころ、それぞれの分割方法に関する特徴的な領域トポロジーが得られたほか、顕著な基底関数依存性や 交換相関汎関数の選択に応じた安定構造の違いによるトポロジーへの反映が見られた。量子エネルギー 密度によって得られた領域境界は運動や力にまつわる量に関する分水嶺として、これらの境界上におけ る局所物性値やその境界面に関する積分量には物理的に大きな意味を含んでいると考えられる。領域分 割によって解析したい主題に合わせて、領域分割方法を選択することが重要であると思われる。 現在、原子軌道展開による周期系について、原子領域や境界を利用して得られる領域エネルギーや領 域双極子モーメント、共有結合次数、局所軌道等を計算するルーチンを開発中であり、領域内や領域間 Li Li 1Li F 1 0.9 (A) 0.8 0.8 0.7 0.7 0.6 0.6 0.5 F Li F 0.4 (A) F 0.2 0.2 Li 1 (A) 0 0.8 (C) F 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 Li Li 0 1 0.7 0.6 0.6 0.5 F (C) (A)(C) F 0.5F Li F F Li 1Li F Li 1 0.9 (A) 0.8 0.7 Li Li (A) (C) 0.9 (A) 0.9 (C) 0.8 0.8 0.7 0.7 F 0.5 (C) 0.3 0.3 0.3 0.2 0.2 0.2 0.2 0.1 0.1 0.1 0 Li 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 Li Li 0 1 0.1 (c) 0.3 0.4 Li 0.5 0.6 0.7 0.8 0.9 Li 1 0.2 0.3 0.4 Li 0.5 0.6 0.7 0.8 0.9 Li 1 0 Li 0 Li (A) Li (A) (C) (A)(C) F F 0.4 0.3 (C) 0.2 (A) (C) F 0.5F Li 0.4 (A) 0.1 0.6 0.6 0.4 0.4 F (b) Li 1Li F F 0.1 Li (a) 0.9 (C) (A)(C) 0.4 (C) 0.3 0 II Li (A) 0.5 0.3 0.1 I Li (A) (C) 0.9 (C) 0.1 F 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 0 Li0 Li 1 0.1 0.2 0.3 0.4 Li 0.5 0.6 0.7 0.8 0.9 Li 1 (d) 図 3 LiF 結晶の原子領域分割(100)および(110) 断面:(a) 基本格子の(100)・(110)面と、最低振動準位 I(3 重縮退)・2番目の振動準位 II(6重縮退)のそれぞれの基準振動モード[12]、(b) 最安定構造での原子領 域、(c) 振動準位 I の基準振動座標上での原子領域、(d) 振動準位 II の基準振動座標上での原子領域。 の相互作用に関する材料系への応用計算をめざしている。原子領域が含有する電子物性解析と振動特性 解析を組み合わせて、大規模動力学シミュレーションへの適用をはかっていく。 参考文献 [1] R. F. W. Bader, Atoms in Molecules: Quantum Theory, (1990), Clarendon Press, Oxford. [2] K. Nakamura, K. Doi, and A. Tachibana, Molecular Regional DFT program package, ver. 1 (2004), Tachibana Lab., Kyoto Univ., Kyoto. [3] K. Doi, K. Nakamura, and A. Tachibana, Periodic Regional DFT program package, ver. 2 (2004), Tachibana Lab., Kyoto Univ., Kyoto. [4] A. Tachibana, “Electronic Energy Density in Chemical Reaction Systems,” J. Chem. Phys., 115, (2001), 3497. [5] J. A. Pople et al., Gaussian 03, (2003), Gaussian Inc., Pittsburgh. [6] F. W. Biegler-König et al., AIMPAC; F. W. Biegler-König, R. F. W. Bader, and T. Tang, “Calculation of the Average Properties of Atoms in Molecules. II,” J. Comput. Chem., 3, (1982), 317. [7] F. W. Biegler-König et al., AIM2000; F. W. Biegler-König, J. Schönbohm, and D. Bayles, “AIM2000,” J. Comput. Chem., 22, (2001), 545. [8] 立花明知・中村康一・坂田健・佐藤貴洋・藤生昌弘・牧野統,領域密度汎関数密度理論計算プログ ラム GREG,京都大学大型計算機センター広報, 34, (2001), 129.(ライブラリプログラム 219/Y7) [9] A. D. Becke, “Density-Functional Thermochemistry. III. The Role of Exact Exchange,” J. Chem. Phys., 98, (1993), 5648. [10] 中村康一・土井謙太郎・中野博史・立花明知,分子系および周期系の領域密度汎関数理論プログラ ム開発,第 9 回理論化学討論会講演要旨集,(2005),p. 124. [11] 中村康一,電子波動関数の摂動論的アプローチによる固体材料の振動特性解析,21 世紀 COE プロ グラム「動的機能機械システムの数理モデルと設計論」平成 16 年度年次報告書,(2005),pp. 159-162. [12] 中村康一,電子波動関数の摂動論的アプローチによる固体材料の振動特性解析,21 世紀 COE プロ グラム「動的機能機械システムの数理モデルと設計論」平成 15 年度年次報告書,(2004),pp. 137-140. Abstract: We have simulated light propagation through a gold thin film with a periodically structured sub-wavelength slit array using a FDTD (Finite Difference Time Domain) method. In the case where the polarization of the incident light is perpendicular to the slit, the light propagation in the slit and the optical near field generation at the slit exit are efficient and the efficiency is almost independent of the film thickness while the efficiency is very small for the case where the polarization is parallel to the slit. The enhancement of the (s-d) electric quadrupole transition of an atom with the optical near field generated at the slit exit is evaluated to examine the possibility of a new optical detection method of the electric quadrupole transition. Key words optical near field, FDTD method, electric quadrupole interaction 1. `³õ-/+"ó½,öp3º0 $î^#%î~`¿0ąî~`% `$wùÁ3íÑ÷cäÛ3qÛ#0*Ąî~` Āô"!ABþ|$` ã·$¯ ´Ä1Ą$}Ê À#øµÂ"ÇÐ"10ą N¢ü±vòªðÍ,ɱmªðÍ"!$` ÌeðÍ%` çðÍ#°'$f¼ß *ĄÝu¶ì 0#+ø2.Ą` ã·#%( 4!À.1"ą Ą` ÌeðÍ%`$z3n£0ü~gñ# .2#ZĄî~`$-#z ` ÆSXÀ0 #-/Ą$f¼ ¦0ąû%1)`$bn# Vî~`3dÀĄl$ü±vòªðÍ$î~`»¥$]`L³¡D9@J#Z 0 3āÃ#è0[1,2] +#Ą$ê×#¾ëÃ"h3P[3]ą %jWýÏ /#î~`3-/f¼Þ¿Ą-/f¼Þü±vòªðÍ3ã· 0¢²3Ö0ą$*Ąt¦Ã";I>@5K63óàÜ#0`$U î~ `$¿# ®_ FDTD ²Finite Difference Time Domain metho#-0¡\åÓ3áą$ Ø©3ÀóàÜî^$ü±vòªðÍȼ$3é[ą 2. 8H;}¨OóàÜ#;I>@Ò÷úQ4t ¦«ï3Ú0ąx1#åÓG?J$NY $¬Ô 3Ëąy¢rĂÕý{Å¢ră#¹ù#õ;I>@ 3Ú0ą;I>@%50 nm#yĄóàÜ$k -&;I>@R÷ú3CHFM= ąa `%³õ3<:7El$s-dü±vòªðÍÙ 3ÿ#685 nm Ą8H;}¨.óàÜ#{Å# a0ą FDTD²%Maxwell¢Î3Ñ÷þ|-&¤÷þ |#chäåÓ¢² 0ą;I>@ 5K6$t¦«ï3â*#x¢r#ø%t¦ Á§T3Àą)Ą®_$åÓ%®_#° 'FGI$ã¸.¥d 0$Ąz¢r$þ|3i cõæą$#-/ĄsoÁ§T3À 1. "UO^*"W¡ =Nd"o«/&P!bj u&"R ®6318"`#~¬\>;! x = z = 5 nm ®ª"s <+! @z/§'GdF"ª]wk/p!` ,)!® 3. Y !!)+q*-¤E"FnmJi"B/® Y>X&-¬\¤,® (A)# x DFGd(B)# y DFGd"]S/®x DF#2934![ DFy DF#l DF,®(a)#2934 (b)(c) #2934V}--250 nm130 nm" ,® 1.2 (A-c) 0.8 300 0.4 200 0.0 100 2 z (nm) (A-b) |E| (a.u) (A-a) 400 0 1.2 (B-c) 0.8 300 0.4 200 0.0 100 2 z (nm) (B-b) |E| (a.u) (B-a) 400 0 0 200 400 x (nm) 600 0 200 400 600 x (nm) 0 200 400 600 x (nm) 2. Transmission Y(B-b)(B-c)!*-,)! y DF"F#2934H/?y IQ2934<¢!#n ª]v- J,®-#y DFGd"]S!#2934H!7:5aZ !cr,®;{x DF"F#e7:5"aZ/N|2934H/?y,®Y(A-b) ")!2934V}`--"2934IQ"ª]#$%+hZ,x] F`¢J ,®Y(A-c)")!¦¨ ,©+S2934*"F£ ,' hZt#f +2934IQ"ª# 1.0 x-pol. F ,® y-pol. 0.8 Y !2934V} 130 nm """ MCat/®"#2934 0.6 *LJ {"Fnm"GdFnm!c, 0.4 ®x DFGd"]S#¤M 0.2 ` ($0g .,® -#") y {T!LJ¥2934! 0.0 , x DFGd"e7:5"aZ/N| 0 100 200 300 400 500 600 ,®&¤M!)+ 20m_K Thickness (nm) ,-#2934H/?y,ª 3. 2934G+QIQ"¦bZ/A,! Transmission Transmission *-§+.-Í_4;=AC4>B*,Ì5@78am i^º¼lX *,ɪ-V-ÍDÌy QRSrº¼Ãt,Ì5@ 78GR<C91\-Í e 4 Á©[ 50 nmÌ5@78w 50 nm 1.0 º¼5@78aOo1ͦ5@ x-pol. 0.8 y-pol. 78i^,̺¼ 0.028 Íx no slit QRSri^º¼5@78G~ÂÄ 200 0.6 nm JEjmÌ º¼ y QRSri^ 0.4 !ɪm Í5@78a 100 nm J 0.2 Eº¼ 90 %1´ÌÁ«©5@781R "#º¼-*-ÍDÌ5@78a 250 0.0 nm ,+m - x QR y QRº¼ 50 100 150 200 250 300 350 ¹¶-Í5@78a 350 nm Ì5@78 Interval (nm) i^*,'º¼s-Í. 4. x QRSri^Ì5@781º¼-RÆi Á«©1º¼-RÆiL¹¶ ,ÌF¨-&§+.-Í y QRSri^Ì5@78at-|̺¼),Z³jm-Í.Ì 5@78DD+0º¼Rµ.-W§+.-ÍÌÁvw5@7 8w_ yÊ5@78a 100 nm qË+s-º¼((jmÌÈh Rº¼<C9%±¸ 1b-Í e 5 5@78U]5@78G~-R 4.0 zySrRzyq-5@78aOo x-pol. 1Í y QRSri^Ì5@78am y-pol. 3.0 /"#DpP-Ì5@78a 2.0 100 nm JE-5@78U]Rzyt m ,Ì'<C9k. 1.0 -nK-Í x QRSrÌe 4 5@78½ 0.0 º¼,Ì5@78U]Rzy'm 50 100 150 200 250 300 350 -5@78aʤ 150 nmË-V Interval (nm) -Í$zy 1 *,'m ̬Ç:?6< 5. B*-ÆijzW*-§+.-Í. *,5@78asi^Ìe 2 ®+.-* Å,^5@78+ÆiÀ,*, ½º¼m DÌ5@78U]Á W* uft1\Rzyt-Í$Ì5@78G~ÂÄ 150 nm * ,m '5@78U]RzytÍ5@78G~ÂÄ 350 nm 5@78U]Rz y *,sÍ *}s5@782A31».-· iR*,ÆcÀn¾ZLM¡ ,¾SrxÇ*-'q.+jm-1°¢Í²¥1Ìm $¯ÆcÀn¾¾RIH0ÆiY¿IHN{ -Í e 6 · iR*-ZLÇ¡,ÆcÀn¾¾ÊPevaËSrxÇi^ÊPincË q-5@78aOo1ÍSrR x QRÍeÌ5@78U]+ z ` 50 nmÌ100 nmÌ200 nm $£dTxg1P1Í*e 5 5@78U] 35 50 nm 100 nm 200 nm eva 25 /P inc 30 20 P 2QP(+)*?Z0JS>]V$k' u(+)*?Z 150 nm YH$ 7 $vu2QP(+)*6<-n3L 2QPK 1hP$ r&"us_@ q\Iogogec 50 nm iB5j 33 1#u2QPlXHFH d$ 7$v %u(+)*6< 6[$sfEH9p'U$ ! $vu(+)*?Z`N.us_@q\ Iogec(+)*6<2QP^Ra` N$v %u(+)*?ZM(+)* 6<sfEq#=9p`N$ mA$vuOC';$iBbu tD,sfE4/K$nWE2T7:$ 8=G$uogecH$v 15 10 5 0 50 100 150 200 250 Interval (nm) 300 350 6. ]PNz_H?y`Jdr z_=\9ZIGCz !&` _H?zy p<z7 + 4.a LVw$6 ]PNz*"ygl/yteu 2 ; FDTD >y0UK Qsr2 wFxx z{ /pmw Gh~yCEy#_HiD,n mwihsr}rz('y]%R`ywMf(8iS 7 150 nm zwjyqz('i4wxmwihsrmz8g|_Hzjx[y 5p]PNzW.y_=\9ZIz(87i 150 nm drv4wxmwihsr mzL8 -^woucCOk_=\9ZI T@pmwzvj]PN:Xz4Y)A 3hypvd M1B [1] S.Tojo, M.Hasuo and T.Fujimoto, Absorption Enhancement of an Electric Quadrupole Transition of Cesium Atoms in an Evanescent Field, Phys. Rev. Lett., 92, (2004) 053001. [2] S.Tojo, T.Fujimoto and M.Hasuo, Precision measurement of the oscillator strength of the cesium 62S1/2b52D5/2 electric quadrupole transition in propagating and evanescent wave fields, Phys. Rev. A, 71, (2005) 12507. [3] S.Tojo and M.Hasuo, Oscillator-strength enhancement of electric-dipole forbidden transitions in evanescent light at total reflection, Phys. Rev. A, 71, (2005) 12508. 微小材料のクリープによる界面はく離発生 工学研究科機械理工学専攻 平方 寛之 Abstract: Crack growth tests are conducted using specimens consisting of a micron- or submicron-thick epoxy film and a glass substrate to clarify the creep crack growth behavior along the interface. The delamination crack grows along the interface due to creep of epoxy and the growth rate is quantitatively evaluated by the in-situ observation. The stress and strain rate fields near the interface crack tip in creep regime are analyzed by finite element method. Using the experimental and analytical results I discus the fracture mechanics parameter which dominates the creep crack growth. In particular, I put the focus on the applicability of the creep J integral which characterizes the intensity of singular stress and strain rate fields ahead of the interface crack tip. Key words: creep, crack growth, delamination, interface, thin film, creep J integral, polymer, epoxy 1. はじめに 低融点金属やポリマーでは,室温下においても熱活性化による原子拡散やクリープが生じるため,これらの 材料は時間に依存する複雑な変形・破壊特性を示す.一方,先進電子デバイスや微小電子機械デバイスは, 寸法がミクロンからサブミクロンメートルオーダーの細線・ドットなどの微小要素から構成されている.これらのデ バイスは,多くの異材界面を有する.界面では変形のミスマッチにより応力が集中するため破壊が生じやすい. しかし,実験的困難さから微小材料の時間に依存する界面破壊特性は明らかになっていない. 薄膜や微小材料の界面破壊は,界面端部からのき裂発生とその伝播によってもたらされる[1-4].発生と伝 播は異なる現象であり,個別に取り扱う必要があるが,本報告ではき裂伝播特性を対象とする.クリープき裂伝 播に関してはマクロな金属材料に関して多くの研究がある[5].これらの材料では,クリープき裂伝播はき裂先 端近傍の応力・ひずみ速度特異場(HRR 場)に支配されていることが明らかとなっており,き裂伝播速度は応 力・ひずみ速度特異場を特徴づけるクリープ J 積分により定量的に評価できる.一方,異材界面においても,き 裂先端近傍に均質材と類似の応力・ひずみ速度特異場が生じることが数値解析により示されている[6].しかし, クリープ界面き裂伝播を対象とした実験に基づく研究は少なく,き裂伝播速度がき裂先端近傍の応力・ひずみ 速度場に律則されているかどうかは不明である.まして,サブミクロン∼ミクロンオーダーの薄膜/基板界面のクリ ープき裂伝播特性は研究されていない.本研究の目的は,サブミクロン薄膜と基板の界面クリープき裂伝播を 支配する力学法則を明らかにすることである.本研究では,その場観察によるき裂伝播特性の観測が可能な 薄膜/基板界面のクリープき裂伝播試験法を開発し,室温でクリープを示す高分子(エポキシ) 薄膜とガラス基 板の界面に適用した.また,エポキシ膜厚を 100 nm ∼ 10 μm オーダーに変えることにより,き裂先端近傍の応 力· ひずみ速度場寸法を変えたクリープ界面き裂伝播試験および力学解析を行い,クリープき裂伝播を支配 する力学量について検討した. 2. 試験方法 図1は,試験片形状および負荷方法を示す[4].高分子材料としては,エポキシ(ペルノックス社製XW-2310) を用いた.本材は,スピンコート法による薄膜化および膜厚の均一化を容易にするため,硬化前に溶剤(メチル エチルケトン)で樹脂の重量比が35%となるまで希釈して粘性を低く設定している(10 mPa s).一方,供試材の基 板にはソーダ石灰ガラスを用いた.本試験片は,硬いガラス基板とシリコン(Si) 基板で高分子薄膜を挟みこむ ことにより,試験中の膜の破壊を防止している.また,ガラス越しに試験中のき裂伝播の様子をその場観察する ことが可能であり,き裂長さを容易に測定できる.厚さ1 mmのガラス基板上にエポキシ溶液を塗布し,スピンコ ートにより製膜した.298 Kで7,200 s風乾させて溶剤を揮発させた.エポキシ膜表面に負荷を伝えるための厚さ 0.3 mmのSi 基板を置き,423 Kのドライオーブン中で3,600 s加熱してエポキシを硬化させた.硬化時の加圧の 有無により膜厚を約0.1 ∼ 0.3 μm (試験片A)と約6 ∼ 8 μm (試験片B)に制御した.試験中における薄膜の破壊を 防止するために,Si基板を接合している部分以外のエポキシを除去した.また,Si基板の破壊を防ぐためにSi 基板背面にステンレス角柱を接着した.試験片を水平に固定し,ガラス端から13 mmの位置(図1における負荷 点)に重量,Ps,の錘を懸垂した.負荷直後から破断に至るまで,光学顕微鏡とデジタルカメラを用いた観察シ ステムにより,ガラス基板越しにき裂伝播の様子を観察した.表1は,それぞれの方法で作製した試験片の寸法 と負荷荷重Psを示す.試験は全て一定温度下(298 K) および大気中で行った. 図 1 試験片および試験システム 表1 No. A1 A2 B1 B2 B3 B4 h, μm 0.1 0.3 6 7 7 8 試験片寸法および負荷荷重 L, mm 1.1 1.2 1.2 1.2 1.2 1.2 W, mm 1.5 1.5 1.8 2.0 2.2 2.0 Ps, N 0.45 0.50 0.80 1.00 1.40 0.90 3. 試験結果 図2は,試験片B4(h = 8 μm,Ps = 0.90 N)の光学顕微鏡写真を経過時間,t,ごとに示す.き裂は,負荷開始 直後に界面端の両角部から発生・伝播し(a),界面端中央部において両角部からのき裂が合体した後(b),試験 片全幅で一様に伝播し(c-d),破断に至った.このき裂伝播挙動は,他の試験でも同様であった.なお,破断面 のオージェ電子分光法による元素分析の結果,き裂はエポキシ/ガラス界面で生じていた.光学顕微鏡写真を 基に,き裂が伝播した面積を試験片幅で除してき裂長さ,a,を求め,a-t関係からき裂伝播速度,da/dt,を計算 した.図3は,da/dtとaの関係を示す.すべての試験において,き裂長さが短いときは,ほぼ一定速度でき裂が 伝播し,その後加速する.膜厚,h,が厚い試験片Aの方がda/dtが速い傾向が見られる. 図 2 界面き裂伝播のその場観察 図3 き裂伝播速度とき裂長さの関係 4. クリープき裂伝播の支配力学量 マクロな金属材では,クリープき裂伝播速度は,定常状態のクリープJ積分J∗によって特徴づけられることが 多くの実験結果により示されている.J∗は,次式によって定義される. ∂u& ⎛ ⎞ J * = ∫ ⎜W *dy − Ti i ds ⎟ Γ ∂x ⎠ ⎝ (1) ε& ij W * = ∫ σ ij dε& ij (2) 0 ここで,W*はひずみ速度エネルギ密度,Tiは積分経路Γ上の表面力, u& i は変位速度,dsは積分経路Γ上の線 素である.また,積分経路Γはき裂先端を取り囲むようにとる.J*は, Γによらず一定の値となり,経路独立性を 有する.しかし,本試験片では,基板による拘束のため,定常状態に達するのに,1014−1016 sのオーダーの時 間を要する.一方,実際の試験時間は104−105 sのオーダーであり,き裂先端はクリープ領域が弾性特異場の 領域に比べて十分に小さい小規模クリープ状態に近い場であった可能性がある.そこでまず,試験結果を基 に有限要素法(FEM)による線形弾性解析を行い,応力拡大係数を求めた.解析には平面ひずみ要素を用い た.解析は,それぞれの試験片の寸法,負荷荷重を用いて試験片ごとに行った.各試験片に対してき裂長さを 変えた複数の解析を行った.用いた弾性定数を表2に示す. 図4 き裂伝播速度と応力拡大係数の関係 表 2 解析に用いた材料定数 Materials Glass Epoxy Silicon Stainless steel E, GPa 70 2.3 130 200 ν 0.20 0.30 0.28 0.30 A, MPa-4.77s-1 2.90 x 10-16 - n 4.77 - 界面き裂先端は混合モード特異応力場であり,総括応力拡大係数Ki = (K12 + K22)1/2 (ここで,K1,K2は複素応 力拡大係数)により特異場の強さを特徴づけることができる.図4は,き裂進展速度da/dtと応力拡大係数Kiの関 係を示す.全ての試験片においてda/dtはKiと良い相関がみられ,Kiの増加に伴いda/dtは増加する.da/dt−Ki 関係は両対数図上でほぼ直線で近似できる.試験片間で比較すると,膜厚のオーダーが等しい試験片A1と A2 (h = 0.1,0.3 μm)および試験片B1-B4 (h = 6 - 8 μm)はそれぞれほぼ同一直線状に分布しており,同じ膜厚 の試験片に対しては負荷荷重によらずda/dtはほぼKiによって特徴づけられる.このことは,本試験の再現性の 高さを示している.しかし,膜厚の異なる試験片間では分布領域が大きく異なる.膜厚の薄い試験片Aは,試験 片Bに比べてKiの変動に対するda/dtの変化が大きい.また,同じKiに対しては膜厚の薄い試験片ほどda/dtの 速い領域に分布している.したがって,Kiによる評価は試験片寸法(膜厚)に依存し,Kiはき裂伝播の支配力学 量として一般性を有していない. 本試験におけるクリープき裂伝播は,定常状態に達する時間より10オーダー以上短い時間で伝播している にもかかわらす,上記のように小規模クリープを仮定した評価も不適当であった.したがって,本試験片のクリ ープき裂伝播は,小規模クリープから大規模クリープに至る途中の遷移状態において生じていると考えられる. 遷移過程であっても界面端近傍ではHRR特異場に類似の応力ひずみ速度場が生じる.また,その特異場の 強さは,遷移過程のクリープき裂伝播を遷移状態でのクリープJ積分,J∗tran,用いて評価可能である.JtranはJ*と 同様に式(1)で定義されるが,その積分経路はクリープひずみが弾性ひずみに比べて支配的な領域内で定義 される.各試験片の遷移状態におけるき裂先端近傍の力学状態をエポキシのクリープを考慮した有限要素法 (FEM)により解析し,J∗tranを求めた.ただし,本研究では,き裂の伝播は考慮せず,各き裂長さのモデルに対し て,そのき裂長さに達するまでの時間が経過したときのクリープJ積分をJ∗tranとした.エポキシのクリープ構成式 はNorton則を用い,材料定数は別途行ったバルク材の単軸クリープ試験で求めた値を用いた.表2に,用いた 材料定数を示す.図4は,da/dtとJ∗tranの関係を示す.傾きは膜厚によって異なっているが,da/dt−J∗tran関係の分 布領域は膜厚に関係なくほぼ一致している.したがって,J∗tranはクリープき裂伝播の支配力学量の有力な候補 であると考えられる.今後,き裂の伝播や一次クリープなどを考慮して遷移過程でのき裂先端近傍の力学状態 を正確に解析することにより,J∗tranの破壊力学パラメータとしての一般性を更に検証する必要がある. 図 5 き裂伝播速度と遷移状態におけるクリープ J 積分の関係 5. まとめ 本研究では,サブミクロン薄膜のクリープによる界面き裂伝播が存在することを示し,その伝播速度を実験的 に評価した.また,き裂伝播速度は遷移状態における界面き裂先端近傍の応力場の強さを表す J∗tran によって 評価できることを明らかにした.クリープによる界面破壊は複雑であり,その総合的理解には,き裂伝播のみな らず,寿命の大きな部分を占めると予測されるき裂発生現象の解明が不可欠である.現在,クリープによるき裂 発生試験法の開発に取組んでおり,今後さらに検討を進める予定である. 参考文献 [1] Kitamura, T., Hirakata, H. and Itsuji, T., Effect of residual stress on delamination from interface edge between nano-films, Engineering Fracture Mechanics, 70 (2003), 2089-2101. [2] Hirakata, H., Kitamura, T. and Yamamoto, Y., Evaluation of Interface Strength of Micro-Dot on Substrate by means of AFM, International Journal of Solids and Structures, 41(2004), 3243-3253. [3] Hirakata, H., Kitazawa, M. and Kitamura, T., Fatigue Crack Growth along Interface between Metal and Ceramics Submicron-thick Films in Inert Environment, Acta Materialia, 54(2006), 89-97. [4] 池田阿希,横山雅憲,平方寛之, 北村隆行,高分子薄膜のクリープによる界面はく離き裂伝ぱ特性,第 55 回理論応用力学講演会講演論文集,(2006),565-566 [5] 平修二,大谷隆一,材料の高温強度論,(1980),オーム社 [6] Biner, S. B., A numerical analysis of time dependent stress fields ahead of stationary interface cracks at creep regime 1 Interface cracks in bimaterials, Engineering Fracture Mechanics, 56 (1997), 513-529. Molecular dynamics simulation of crack initiation at bi-material interface edges 工学研究科機械理工学専攻 Fulin Shang Abstract: This study examines the critical conditions of crack initiation from a range of interface edges with different angles. Molecular dynamics simulations are performed to check the onset of fracture at the interface edges with angles 45°, 90° and 135°, respectively. As a first step, the case of perfectly brittle fracture is addressed. It is found that: (a) the maximum stresses reach the ideal strength of the interface for the three models, at the onset of fracture; (b) the fracture interface energies are prescribed by the cohesive interface energy, which is an intrinsic property of the interface material. Thus, unified criterion for the crack initiation at these edges could be formulated through either the ideal strength or the cohesive energy of the interface. And the crack initiation at the interface edges is actually controlled by the maximum stresses or the cohesive interface energy. Key words: crack initiation, interface edge, molecular dynamics simulation, stress, cohesive energy 1. Introduction This work continues the study begun in H.16 and also supported by this COE Program, where the correlations between fracture criteria for crack initiation at two interface edges and that for crack propagation along an interface crack were investigated in particular. The objective of this study is to examine crack initiation from a range of interface edge with different angles, and see if a unified criterion for the crack initiation at these edges could be formulated. Further, this study will focus on the interface edges between thin film and thick substrate, where films are typically of sub-micron meter thickness or even thinner. The first issue to be considered is perfectly brittle fracture of the interface. For this purpose, molecular dynamics simulations are conducted to check the conditions of the onset of fracture at the interface edges with a variety of angles. 2. Simulation method and models In principle, the same simulation method as in the H.16 study is used here. That is to say, Morse pair potential describes the interaction between atoms, Φ ( rij ) = D[e −2α ( rij − r0 ) − 2e −α ( rij − r0 ) ], (1) where Φ is a pair interaction between atoms i and j whose separation is given by rij, D, α and r0 are the potential parameters. A combination rule for the interactions between two dissimilar materials 1 and 2 is adopted to introduce a weak interface between materials, as follows, D12 = 101 ( D1 D2 )1 2 , α12 = 12 (α1 + α 2 ) , r012 = (λ1λ2 )1 2 + ln 2 α 12 , λi = r0i − ln 2 α i , i = 1,2 . (2) More details about the simulation method please refer to the H.16 Research Report [1]. In this study, three sets of simulation models are considered, see Figure 1: (a) interface edge with angle 45°; (b) interface edge with angle 90°; (c) interface edge with angle 135°. As in the H.16 study, the outer boundaries of these models are displaced stepwise with an increment du according to the linear elastic FEM displacement fields. The MD cells of the three interface edges are stretched until the onset of fracture occurs at the interface edges. Then the stress fields and cohesive energy along interface are examined in detail. 3. Simulation results and discussions 3.1. Stress fields Figure 2 shows the stress distributions along the interface of the model (a) during the loading steps, where both normal and shear stresses are collected. It is clear that the shear stresses are quite low comparing to the normal ones, and the cell is mainly under tension deformation. Figure 3 gives a close-up view of the normal stress results. By observing the stress plot at the apex of the interface edge, we could see that the normal stress starts decreasing at the step of du=7. From the definition of fracture proposed in the H.16 Report, therefore, one conclusion could be reached as follows, i.e. the crack initiation occurs at the step of du=7. The snapshots of atoms for this model at this critical step and after this step are shown in Fig. 4. load load Material 1 10 μm Material 1 10 μm 36 μm 36 μm interface interface 36 μm 36 μm Material 2 10 μm Material 2 10.8 nm MD region 10.8 nm MD region 10 μm 10.8 nm Material 1 10.8 nm Material 1 Material 2 z Material 2 z x y periodic x y periodic 18 nm 36 μ m 18 nm 36 μ m (a) Interface edge with angle 45o (b) Interface edge with angle 90o load Material 1 10 μm 36 μm 10.8 nm Material 1 interface 36 μm 10.8 nm MD region 10 μm Material 2 Material 2 z x y periodic 15.2 nm 36 μ m (c) Interface edge with angle 135o Fig. 1. Simulation models For the model (b), the stress results along the interface during the loading steps are shown in Fig. 5, including both normal and shear stresses. Similarly, the shear stresses are relatively much lower than the normal ones, and the cell is also mainly under tension deformation. Figure 6 gives a close-up view of the normal stress results. We could observe that the normal stress at the very apex of the interface edge starts decreasing at the step of du=10. Therefore, the conclusion is that the crack initiates at the step of du=10 for the model (b). The snapshots of atoms for this model at this critical step and after this step are shown in Fig. 7. For the model (c), the stress results along the interface during the loading steps are displayed in figure 8, including both normal and shear stresses. Similarly, the shear stresses are much lower than the normal ones, and the cell is again mainly under tension deformation. Figure 9 is a close-up view of the normal stress results. We could observe that the normal stress at the very apex of the interface edge starts decreasing at the step of du=19.5. Therefore, the conclusion is that the crack initiates at the step of du=19.5 for the model (c). The snapshots of atoms for this model at this critical step and after this step are shown in Fig. 10. 6 6 du=7 du=8 5 du=7 du=9 du=6 du=6 4 du=5 3 σzz du=4 σzz , GPa Stress, GPa 4 2 du=2 1 du=5 du=4 3 2 1 du=2 σxz 0 du=8 5 0 du=0 −1 −1 6 8 10 12 14 16 18 6 8 10 12 14 16 18 Distance from the apex of interface edge, nm Distance from the apex of interface edge, nm Fig. 2. Stress distributions along the interface for edge with angle 45o Fig. 3. Normal stress distributions along the interface for edge with angle 45o (a) du = 7 (b) du = 9 Fig. 4. Snapshots of atoms for the interface edge with angle 45o 6 6 du=9 5 du=10 du=11 3 du=12 σ zz du=8 σ zz, GPa Stress, GPa 4 du=7 2 1 4 du=8 du=7 3 du=6 du=5 2 du=4 du=2 0 0 −1 −1 8 du=12.25 1 σ xz 6 du=10 du=11 du=12 du=9 5 10 12 14 16 18 6 8 10 12 14 16 18 Distance, nm Distance, nm Fig. 5. Stress distributions along the interface for edge with angle 90o Fig. 6. Normal stress distributions along the interface for edge with angle 90o (a) du = 10 (b) du = 12 Fig. 7. Snapshots of atoms for edge with angle 90o With these results in hand, we can now compare them with the maximum stress that the interface material can afford, that is, the ideal strength of the interface. Figure 11 shows this comparison graphically, where the solid line denotes the ideal strength. It is seen that the stress values for the three models studied here fall approximately around the curve of the ideal strength. The first two edges fractured at the stress levels very close the ideal strength, while the deviation from the ideal strength for the interface edge model (c) is relatively larger. But the three interface edges all tend to fracture when their maximum stresses rise to the ideal strength. This implies that the onsets of the fracture at the three interface edges are controlled by the ideal strength of the material. This conclusion is understandable, if one recalls the results in H.16 study. 3.2. Cohesive energy Figures 12, 13 and 14 collected the simulation results of the cohesive energies Ei along the interfaces of the three models, at the onset of fracture determined according to the above stress calculations. Also, the cohesive energy of the interface material Ei,cohesive is inserted into these figures, which is represented by the dashed line. It is understood that this cohesive energy is an intrinsic property of the interface material itself. From this comparison, we could see that: (a) the interface energies at the very apexes of the three edges are just lower than the cohesive energy, at the onset moment of fracture; (b) while the other portion of the interfaces are still held with the energies that are higher than the cohesive energy. The critical condition for crack initiation at the interface edges with different angles, therefore, can also be unified through the cohesive energy. That is, Ei (interface edge) ≤ Ei,cohesive . (3) 4. Conclusion This study examined the conditions of crack initiation from a range of interface edges with different angles. Molecular dynamics simulations were performed to check the onset of fracture at the interface edges with angles 45°, 90° and 135°, respectively. As a first step, the case of perfectly brittle fracture was addressed. It is found that: (a) the maximum stresses reach the ideal strength of the interface for the three models, at the onset of fracture; (b) the fracture interface energies are prescribed by the cohesive interface energy, which is an intrinsic property of the interface material. Thus, unified criterion for the crack initiation at these edges could be formulated through either the ideal strength or the cohesive energy of the interface. And the crack initiation at the interface edges is actually controlled by the maximum stresses or the cohesive interface energy. References [1] Shang, F., Correlation between fracture criteria for interface crack initiation and propagation: an atomistic simulation, H.16 COE Research Report, Kyoto University, 2004. 4.5 4.5 4 4 3.5 3.5 3 Normal Stress, GPa du=21.5 2.5 du=19.5 normal stresses 2 1.5 du=16.5 du=12 1 0.5 0 −0.5 6 8 10 du=21.5 du=19.5 2 1.5 du=16.5 du=12 1 0.5 shear stresses 4 3 2.5 12 14 0 16 4 6 8 10 12 14 16 Distance, nm Fig. 8. Stress distributions along the interface for edge with angle 135o Fig. 9. Normal stress distributions along the interface for edge with angle 135o (a) du = 19.5 (b) du = 21.5 Fig. 10. Snapshots of atoms for edge with angle 135o 3.5 10 3 9 2 ideal strength Interface energy, J/m Shear stress, GPa 2.5 2 1.5 1 edge 90 edge 45 0.5 du= 2 4 5 6 7 8 8 9 7 0 1 2 3 4 11 5 4 5 3 6 2 6 12 edge 135 0 6.57 J/m 10 7 8 9 10 11 Normal stress, GPa 12 13 14 15 16 17 18 Distance, nm Fig. 11. Comparison of the maximum stresses for three models with the ideal strength Fig. 13. Interface energies for edge with angle 90o 9.5 10 du=2 9 4 5 8 7 6 6 7 6.57 J/m Interface energy, J/m 2 Interface energy, J/m 2 9 2 5 8 4 du=9 3 8 7.5 7 6.57 J/m 2 6.5 2 1 8.5 6 7 8 9 10 11 12 13 14 15 16 17 18 Distance, nm Fig. 12. Interface energies for edge with angle 45o 4 6 8 10 12 14 16 Distance, nm Fig. 14. Interface energies for edge with angle 135o MIT 武者修行成果報告 工学研究科機械理工学専攻 久島 祥嘉 1. はじめに 私は,21世紀 COE プログラム「動的機能機械システムの数理モデルと設計論」-複雑系の科学による機械 工学の新たな展開―における若手武者修行制度の助成を受け,2006年1月23日から2月17日までの約4週 間アメリカ合衆国マサチューセッツ工科大学(MIT)に滞在した.今回の滞在目的は同大学の Sydney Yip 教授 のグループを訪問し,材料のマルチスケールモデリングに関する最新の情報収集および,議論を交わすこと であった.Yip 教授は材料の物理特性に関して,電子・原子レベルから連続体までの広いサイズスケールから 解明するという研究を行っている. 本報告では MIT 滞在中に行った活動についてまとめる. 2. マサチューセッツ工科大学滞在 2.1 マサチューセッツ州 マサチューセッツ州は600万人以上と全米代13位の人口を持ち,「大きな丘の場所」を意味する先住民の言 葉「マサチューセッツ」から名づけられた.これはボストンの南,ミルトンにある小さな山にちなんだものである. 1620 年にメイフラワー号でピルグリムファーザーズがプリマスに入植し,その後に清教徒達が続き,マサチュー セッツ湾殖民地を建設した.清教徒たちは信仰の自由を求めてマサチューセッツに来たものの、彼らは自分た ちのもの以外の宗教には寛容ではなかったため、アン・ハチンソン、ロ ジャー・ウィリアムズ、トマス・フッカーの ような人々はマサチューセッツを離れ、南に向かった。ウィリアムズはロードアイランド植民地を建て、Hooker は コネチカット植民地を建てた。アメリカ独立戦争の初期の戦いのいくつかはマサチューセッツで起こり,その中 にはレキシントン・コンコードの戦い,バンカーヒルの戦い,ボストン包囲戦も含まれる.そのため多くの史跡が 残っている.清教徒の影響が大きいマサチューセッツ州では州法により、日曜日アルコール類を店で販売する ことは禁止されている. 2.2 ボストン&ケンブリッジ ボストン近郊には60以上の大学があり平均年齢は 26.5歳で ある.東海岸に位置し,ニューヨークにも近く比較的安全な 街として知られている.有名大学が数多くあるため世界中か ら多くの学生,研究員が訪れている.学生が多く在住してい るためアメリカの都市のわりには比較的安全な場所として知 られ,快適な生活を送ることができる.また,独立戦争の舞 台となったところでもあり多くの史跡を見ることもできる.ボス トンを本拠地として活動している世界でもっとも有名なオー ケストラのひとつであるボストン交響楽団がある.港町ボスト ンを象徴するウォーターフロントは空路の発達とともに閑散 とした倉庫街となったが,1960 年代からの再開発計画により 図1:チャールズ川とボストンの夜景 観光名所としての賑わいを取り戻している.ここから様々な クルーズが出航しており特に夏期のホエールウォッチングが人気である.またここには 1773 にイギリス政府の 課税に反対した市民が茶箱を投げ捨てた有名なボストン茶会事件が起こった3隻の船のうちの1隻の複製船が 停泊している. ボストンの北側のチャールズ川(図1)を挟んで対岸に位置するケンブリッジは学園都市として名高くハーバード 大学やマサチューセッツ工科大学などがある.ハーバード大学付近は多くの劇場,カフェ,レストランが集まっ ており夜遅くまで人通りが多い.一方,MIT 周辺はハーバード大学とは対照的に近代的なビルに囲まれてい る. 2.3 マサチューセッツ工科大学 マサチューセッツ工科大学はボストンの北側を流れるチャール ズ川の北側の学術都市ケンブリッジにあり MIT の略称で知られ ている.MIT は世界のハイテク技術をリードする工科部門では屈 指の大学である.同大学の博物館では小規模だが MIT の研究 成果を知ることができる.ホログラフィー,機械工学,人間の目で は見ることのできないスローモーションなどの展示のほか,“MIT いたずらの殿堂(MIT Hall of Hack)”という今まで MIT の学生が 図2:金メダルがかけられたロジャースビル 実際に行ったいたずらの数々が展示されている.イタズラの中に は非常に高度なものもあり研究以外のことにも頭をつかって楽しもうという MIT の学生の姿を垣間見た気がし た.このように日常生活でも遊び心を忘れないところが優れた研究につながっているのではないかと感じた.ト リノオリンピックの金メダルがかけられた MIT の象徴,ロジャースビル(図2).MIT の建物はすべて地下通路で つながっており,一度も外に出なくても建物間を移動できる.滞在先研究室の学生に地下通路を案内してもら ったが,地下では外の風景がまったくわからないため自分がどこにいるのか迷ってしまうそうだ.地下通路とい っても暗く細いものではなく建物内の廊下と同じような感じで地下に研究室もある.通路には迷わないように地 下通路の地図が様々なところに貼り付けられていた.今回の滞在では Department of Nuclear Science & Engineering の Sydney Yip 教授のグループを訪問した.同グループでは耐火ホウ化物の高温酸化,分子の電 気伝導度,分子アクチュエータなど幅広い分野の研究が行われている.また,筆者の滞在期間中に Yip 教授 は 70 歳の誕生日をむかえられ,研究室のメンバー達とともにささやかなお祝いをした.研究室の一席を借りる ことができ,同室のポスドクの Clemens, Xi 博士と研究の分野だけでなく互いの国の文化など様々なことについ ても意見を交わすことができた.滞在5日目には研究発表の場をいただき,これまでの私の研究成果を発表し た.数多くの質問,意見をいただきまた多くの議論もでき予定時間を大幅にうわまわり,研究員・学生の豊富な 好奇心,基礎知識に驚き,学ぶところが多くあった.ポスドクの Xi との議論の中から新たな研究テーマができ, 滞在中に共同研究を行うことができた.以下にその概要を示す. 3. Si/Ge 界面 Si/Ge 界面を有する構造体では,格子定数の差から 生じるひずみによるバンド構造の変化および Si,Ge の Ge 価電子バンドのオフセットにより生じる高いホール移動 度など優れた電気的特性を示し,ナノデバイスの要素 材料として注目されている.また,近年製造が可能にな Si った Si/Ge 界面構造を有するナノワイヤでは非常に高 い電気伝導度が観察されている[1].これまでに Si/Ge 界面構造の電気的特性に関する多くの解析がされてき ているが[2-5],それらの解析では Si/Ge 界面において 格子定数のミスマッチにより等方ひずみが生じるとして 5 nm その価電子バンドオフセットやバンドギャップを計算し ている.Si 基板上に Ge 薄膜が堆積しているような2次 図3:Si/Ge ナノワイヤの断面図[1] 元構造においては Si, Ge に加わるひずみは界面上で 等方的であるが,Si/Geナノワイヤ等のような1次元構造 の場合では界面には非等方ひずみが生じることが考えられる.本研究では Si/Ge 界面上の非等方ひずみが電 気的特性に与える影響を第一原理計算を用いて厳密に評価した. 3.解析モデル 図2に Si/Ge 界面構造の解析モデルを示す.(100)界面を有している.ユ ニットセルは Si,Ge 原子を8個ずつ含み x,y,z 方向に周期境界条件を適用 した.解析には第一原理計算ソフトの VASP[6,7]を使用し,擬ポテンシャルは PAW[8],交換相関項には一般化密度勾配近似(GGA)[9]を用いた.カットオ フエネルギーは400eV とした.x,y 方向の格子定数を Si から Ge までの間で 様々な組み合わせで変化させ原子構造および z 方向の応力が0になるよう に z 方向のセルサイズを緩和し,バンドギャップエネルギーを計算した.緩和 計算には共役勾配法を用いた. 4.解析結果 図4に x,y 方向の格子定数 ax,ay とバンドギャップエネルギーEgap の変化を 示す格子乗数が Si から Ge に変化するにつれバンドギャップが減少するが, ax=ay の時,すなわち等方ひずみの場合はバンドギャップの減少は少なく 図4:解析モデル 格子定数が Si から Ge に変化してもバンドギャップは存在し半導体の性質を 示す.一方,ax=0 または ay=0 のときは Egap の減少の度合いが大きく Si から Ge の格子定数の範囲内でバンドギャップが 0 になり半導体から導体に遷移したことがわかる.今回の滞在中に 行った解析はここまでであるが,今後もこの課題について Xi と議論をしながら研究を進め,4月中には論文に まとめて投稿する予定である.そのために,なぜ非等方ひずみに対してバンドギャップの減少が大きいのかを 理論的に説明する必要があると考えられる.まずひずみが加えられた際の原子配置の変化を詳細に観察し, 等方ひずみの場合と比較する.また,Si/Ge 界面の重要な性質として,Si の価電子バンドの最大値と Ge の価電 子バンドの最大値との差,価電子バンドオフセットがある.この差が Ge 内に正孔をつくり電気伝導度を増加さ せる原因になっているのではないかという議論がなされている.そこで,バンドギャップの変化だけでなく価電 子バンドオフセットに対しても解析を行い,Si/Ge 界面構造の電気的特性に対してより深い考察を行う. ay 半導体 a Ge a Si a Si a Ge a Ge a Si 導体 図5:非等方ひずみに対する Si/Ge 界面のバンドギャップエネルギー変化 a Ge ax 5.最後に 今回の滞在中は,MIT の Yip 教授のグループを訪問し,これまでに私の行った研究および最新の研究テーマ に関する建設的な議論をすることができ,また共同研究を行うことができ成果を残すことができるなど非常に有 益なものであった.またポスドクの研究員,同じ博士課程の学生と交流をもつことができ意見を交換することが できたことも大きな収穫であった.彼らと知り合えたことで今後も積極的に研究に対する意見交換をすることで 将来の自分の研究の発展につなげられると考える.学生や研究員達の知的好奇心の大きさ,疑問を持ったら どんなに小さなことでも考え,質問する積極性に非常に感銘し,学ぶところが大きかった.今回の MIT での滞 在は,自身の視野を広げ,研究に役立つ多くの情報をえることができ非常に有意義なものであった.今回得た ことを今後の研究活動に生かしていきたいと考えている. 6.謝辞 今回,若手武者修行制度により筆者に海外での貴重な経験を与えてくださいました,21世紀 COE プログラム 「動的機能機械システムの数理モデルと設計論」-複雑系の科学による機械工学の新たな展開―拠点リーダ ー航空宇宙工学専攻 土屋和雄教授,材料グループリーダーであり,筆者の指導教官でもある機械理工学専 攻 北村隆行教授をはじめとする関係者の皆様に厚く御礼申し上げます.今回の滞在で筆者を受け入れてく ださいましたマサチューセッツ工科大学 Sydney Yip 教授および研究室の皆様,同グループで共同研究を行う に際して有益なアドバイスをいただきました Xi、Clemens, Xeofeng に心より感謝いたします.また,筆者が滞在 先から計算機にアクセスできるように設定をしてくれた筆者所属研究室の修士課程の学生 嶋田隆広に感謝し ます. 参考文献 [1] W. Lu, J. Xiang, B.P. Timko, Y. Wu, and M. Lieber, One-dimensional hole gas in germanium / silicon nanowire heterostructures, PNAS, 102 (2005), 10046. [2] C. G..Van de Walle, Theoretical calculations of heterojunction discontinuities in the Si/Ge system, Phys. Rev., B34 (1986), 5621. [3] M. Ikeda, K. Terakura, and T. Oguchi, Theoretical study on the electronic structure of (Si)m/(Ge)n superlattices, Phys. Rev., B48 (1993), 1571. [4] S. Satpathy, R. M. Martin, and C.G.Van de Walle, Electronic properties of the (100) (Si)/(Ge) strained-layer superlattice, Phys. Rev., B38 (1988), 13237. [5] L. Colombo, R. Resta, and S. Baroni, Valence-band offsets at strained Si/Ge interfaces, Phys. Rev., B44 (1991), 5572. [6] G. Kresse and J. Haffner, Ab initio molecular dynamics for liquid metals, Phys. Rev., B47 (1993), 558. [7] G. Kresse and J. Furthmüller, Efficient interactive schemes for ab initio total-energy calculations using a plane-wave basis set, Phys. Rev., B54 (1996), 11169. [8] G. Kresse and J. Joubert, From ultrasoft pseudopotentials to the projector augmented wave method, Phys. Rev., B59 (1999), 1758. [9] J. Perdew and Y. Wang, Accurate and simple analytic representation of the electron-gas correlation energy, Phys. Rev., B45 (1992), 13244. 第一原理計算に関する武者修業成果報告 工学研究科機械理工学専攻 1. 木下 佑介 はじめに 私は,21 世紀 COE プログラム「動的機能機械システムの数理モデルと設計論」―複雑系の科学 による機械工学の新たな展開― における若手武者修業制度の助成を受け,2006 年 2 月 26 日から 3 月 9 日までの 12 日間ヨーロッパに滞在した.今回の滞在目的は,第一原理計算による材料の変 形特性評価に関する研究および第一原理計算アルゴ リズムとそのコーデ ィングについての最新の 研究成果に触れ,議論を交わすことであった.まず,チェコ・ブルノのマサリック大学の Mojmir Šob 教授の研究グループを訪問した.Šob 教授は,第一原理計算を用いて,Fe や Ni といった磁性 材料の強度や変形に伴う磁性変化に関しての研究を行われている.また,ブルノ滞在中には,Šob 教授の下で学ぶ学生の博士論文公聴会がオーストリア・ウィーンで行われ,それに同行した.次 に,ベルギー・ルーヴァン -ラ- ヌーヴのルーヴァンカトリック大学の Xavier Gonze 教授の研究グ ループを訪問した.Gonze 教授は第一原理計算ソフト ABINIT[1]-[3] の開発リーダーであり,日々 その開発に精力的に取り組まれている.筆者も ABINIT を使用して研究を行っている.本報告で は,今回のヨーロッパ滞在での活動およびそれから得たこと・感じたことを記す. 2. マサリック大学滞在 2.1 ブルノ ブルノはチェコの南東方に位置し,チェコの首 都プラハから直線距離にして約 180km である.人口は約 38 万人で,チェコ第二の都市である.ブルノは三つの方面から 山に囲まれ,南方だけには平野が広がる.その平野にはスヴ ラットカという旧市街からちょっと離れている川が流れてい る.町自体には様々な丘があり,ブルノは要塞を創設するた めに完璧な場所だったと言える.旧市街の構造は他の古い町 と比べて珍し く,昔からの外郭がちゃんとあるが,その真中 図 1: シュピルベルク城 には一つの広場ではなく沢山の広場がある.ブルノ中央駅前 には大きいトラム停があり,トラムで色々な名所まで行ける. 筆者も移動にはトラムを頻繁に利用した.駅からマサリコヴァ 通りを北へ 5 分程歩くと,自由広場という一番範囲が広い広 場に着く.ブルノのほとんど 全ての名所はその周辺に密集し ている.この町の誇りであるシュピルベルク城 (図 1) は旧市 街の西方の丘の上にある.どこからも見える城と二つの大き い塔のある聖ペトロ聖パウロ教会 (図 2) のおかげで町を迷わ ず歩くことができる. チェコはビールで有名であり,国民一人当りのビール消費 量は世界一である.ビールの種類も豊富で,街ごとにあるい は店ごとに異なるビールが存在し,ブルノにも”Starobrno” と いうビールがある.ブルノ滞在中は,いろいろなビールを味わ うことができた (図 3).驚いたのが多くの人が昼間からビール 図 2: 聖ぺテロ聖パウロ教会 をたくさん飲むこととその安さである.種類にもよるが,100 円程度で 0.5l のビールを飲むことができ,しかも水よりも安いことには驚愕した. 2.2 マサリック大学 マサリック大学 (図 4) は,いくつかの学部・学科がブルノ市内に点在している が,いずれも街の中心部からはそれほど離れておらず,徒歩またはトラムで 10 分もあれば到着する. 今回は,Department of Theoretical and Physical Chemistry の Mojmir Šob 教授の研究グループを訪問した.Šob 教 授は,第一原理計算による材料の理想強度解析に関する多く の研究論文を発表され,その他にも,第一原理計算と CALPHAD 法を組み合わせた二元系・三元系状態図の理論計算や, 相変態に伴う Fe3 Al,Ni3 Al など の磁性変化に関する第一原 理解析といった研究を行われており,その研究範囲は多岐に 渡る.Šob 教授は,現在,学科長を務めておられることから, 図 3: 乾杯する筆者と Šob 教授 Šob 教授の研究室だけでなく,学科内のほとんど 全ての研究 室を訪問し ,スタッフの方々に今どのような研究を行ってい るのかを丁寧に説明していただくことができたことは非常に 幸運であった. 滞在中は,研究室の一席を借りることができ,同室の博士 課程の学生 Martin と Tomáš とは研究のことだけでなく,互 いの国の違い等についても意見を交わすことができた.特に, Tomáš は,日本についてとても興味があるようで, 「 こんな時, 図 4: マサリック大学 日本語では何と言うのか」などと度々質問された.このよう に互いの文化に触れ,交流を深めることができたことも今回 の滞在で得た大きな収穫の一つであった. 滞在 3 日目には,研究発表の場を設けていただき,これまでの研究および現在行っている研究 について約 40 分程発表した.その後の質疑応答では,こちらが想定していた質問もあったが,そ れ以上に思いもよらない質問が多く,自身の研究を見つめ直すのに非常に有益な場であった.終 了後,個人的に質問に来る方も何人かおり,自身の研究に興味を持ってくれていることを非常に うれしく感じた. 3. ウィーン大学訪問 3.1 ウィーン ウィーンはオーストリアの都市名および州名で あり,同国の首都でもある.人口は約 163 万人で,クラシック音 楽が盛んなことから「音楽の都」 ・ 「楽都」などと呼ばれる.現在, 日本人指揮者小澤征爾がウィーン国立歌劇場の音楽監督を勤めて おり,おそらく彼は現在,ウィーンで最もよく知られている日本 人の一人である.事実,ウィーンでお会いした方のほとんどは,筆 者が日本人であることが分かると,”Do you know Seiji Ozawa?” 図 5: ボルツマン通り といった感じに質問してきた. ウィーンの滞在時間は実質半日しかなかったため,市内を見物す る時間はほとんどなかったが,移動に用いたトラムの中からウィー ン国立歌劇場や市庁舎などいろいろな歴史的建造物を眺めること ができ,窓越しではあったがその圧倒的存在感を感じることがで きた. ウィーンが音楽の都であること,音楽に留まらず芸術が盛んであ ることは日本でも良く知られているが,一方でウィーンは偉大な研 究者を多く輩出した都市でもある.徒歩で移動中,”Boltzmann” という単語が目に入ったので眺めていると,Šob 教授が, 「 この通 りはボルツマン通りといってあの有名な物理学者の名前に由来し ているのだよ」と説明して下さった (図 5).他にも精神分析で有 名なフロイト (図 6) など ,ウィーンに縁のある研究者は多いとの 図 6: フロイトが暮らした場所 ことであった. 3.2 博士論文公聴会 ウィーンを訪れた目的は,Šob 教授の学生がウィーン大学で行う博士論文 公聴会に参加することであった.講演は,第一原理計算を用いた固体材料の機械的特性のモデル 化と題し ,具体的には原子レベルからき裂および転位を対象とする研究であった.き裂および転 位をいくつかの原子からなる単純なモデルにより模擬し ,進展抵抗に相当する新たに定義した物 理パラメータを第一原理計算により求めると,その値は材料に依らずほぼ一定となるというもの であった.材料の破壊を体系的に取り扱う学問として破壊力学というものが存在するが,この研 究は破壊現象を原子レベルから体系化しようとするものであると言える. 公聴会は,講演が約 40 分,質疑応答が約 30 分で,質疑応答終了後は主査・副査のみ教室に残り 約 15 分程度の審議が行われるというものであったが,その流れや雰囲気は筆者がこれまでに参加 した京大でのそれとほとんど 同じであった.審議中に講演者に公聴会の感想を尋ねたところ, 「最 初の質問が一番厳しかった」「 ,これを質問されたらど うしようかと思った」などの答えが返ってき たが,彼は若干興奮気味で,その顔は緊張から開放されたせいか非常ににこやかであった.公聴 会終了後の開放感のようなものが万国共通であることを感じた瞬間であった. 4. ルーヴァンカト リック大学滞在 4.1 ルーヴァン -ラ-ヌーヴ ルーヴァン -ラ-ヌーヴ (図 7) は, ベルギーの首都ブ リュッセルから約 23km 東方に位置する学 園都市であり,その歴史は 50 年程度である.都市というより も大学を中心とした集落と言った方が適しているかもしれな い.駅前広場を中心に,レストランや食品雑貨店,生活用品 店が並び,不自由なく生活することができる.ベルギーの名 物であるワッフル (図 8) を販売する露店もあり,筆者も食べ てみたが非常に美味であった.ベルギーはブ リュッセル辺り を境に南北で使用言語が異なる.北側はオランダ語の方言の 図 7: ルーヴァン -ラ- ヌーヴの駅前 一種であるフラマン語,南側ではフランス語が使用されてお り,ルーヴァン -ラ- ヌーブはフランス語圏である.ベルギーに はルーヴァン (ルーヴェン ) という都市もあり,むしろこちら の方が有名である.筆者はホテルの予約等をするまでは,ど ちらも同じ都市のことを指していると勘違いしていたのだが, 旅行会社の方の指摘で全く別の都市であることを知った.こ のことを滞在の世話をしてくれた博士課程の学生 Simon に話 したところ,やはり勘違いしてルーヴァンに行ってしまうゲス 図 8: ベルギー名物ワッフル トも多く,度々ルーヴァンまで迎えに行くとのことであった. 4.2 ルーヴァンカト リック大学 (Université Catholique de Louvain ; UCL) ルーヴァンカトリック大学 (以下 UCL) は,ルーヴァン -ラ-ヌーヴ内にあ るというよりも,UCL が街を形成していると言っても過言ではない.そのため,どこまでが大学 構内でどこからが街なのかがはっきりせず,慣れていないとかなり迷ってしまう.筆者も初日は 中々目的地に辿り着くことができなかったが,親切な女性が地図を用いて丁寧に説明して下さっ たおかげで何とか辿り着くことができた. UCL 滞在の主目的は,Unit of Physico-Chemistry and Physics of Materials (PCPM) の Xavier Gonze 教授の研究グループを訪問し ,Gonze 教授が開発リーダ ーを務めておられる第一原理計 算ソフト ABINIT に関して議論することであった.Gonze 教授はお忙しい中,筆者の質問一つ 一つに丁寧に答えてくださり,現在の開発状況や開発の目指す先など についても詳し く説明し てくださった.中でも筆者の興味を引いたのが,GW 近似計算の更なる充実とウェーブレット変 換導入による計算の高速化である.第一原理計算の根幹を成す密度汎関数理論 (Density Functional Theory ; DFT) はバンドギャップを実験値の半分程度に過小評価してしまうことが知られ ており,これを解決する手法が模索されている.GW 近似もそのような手法の一つであり,GW 近似計算を備えていることが ABINIT の大きな強みの一つである.これを更に拡張することに より研究の幅が拡がることは間違いない.第一原理計算では,実空間と逆空間とのやり取りを 行うが (実空間法を除く),これまではフーリエ変換によってそれを行ってきた.一方,ウェー ブレット 変換も周波数解析法の一種であり,全デ ータ量を N とするとその計算量は O(N ) で ある.これは計算量が O(N log N ) であるフーリエ変換よりも高速であることを意味し ,これ を導入すれば 計算の高速化が図れる.計算高速化によってもたらされる恩恵は言うまでもない. PCPM には,Gonze 教授率いる ABINIT 開発グループの他 に二つのグループがある.一つは Jean-Christophe Charlier 教 授が率いるカーボンナノチューブの電気的特性に関する第一 原理解析を行うグループ,もう一つは Gian-Marco Rignanese 教授率いる High-k 材料に関する第一原理解析を行うグループ である.筆者も主に機械的特性ではあるがカーボンナノチュー ブの研究を行っているので,Charlier 教授との議論はかなり 有意義なものであった.Rignanese 教授は,知識がなくても 分かるようにゆっくりと丁寧に研究内容を説明してくださり, 図 9: PCPM の組織図 PCPM の組織図 (図 9) まで描いてくださった.この図を用い て誰がどこに所属していてどんなことをしているのかをあら かじめ説明していただいたことが,後に各スタッフ・学生と議論する際,非常に役に立った. UCL でも研究発表の機会を与えられ,マサリック大学で行ったものとほぼ同じ 発表を行った. 一度話したことがあるせいか,以前よりも緊張せずに発表を行うことができた.ここでも多くの 有用な質問・コメントをいただくことができ,今後の研究に多いに役立つと思われる. 5. おわりに 今回,12 日間という短い期間ではあったが,マサリック大学,ウィーン大学,ルーヴァンカト リック大学を訪問し, 「 第一原理計算による材料の変形特性評価に関する研究および第一原理計算 アルゴ リズムとそのコーデ ィングについての最新の研究成果に触れる」という目的を達成するこ とができた.また,各研究室のメンバー,特に同じ博士課程の学生と交流できたことも大きな収 穫であった.今回のヨーロッパ滞在は,自身の視野を広げ,研究に役立つ多くの情報を得ること ができたいう点で非常に有意義なものであった.今回,得たこと・感じたことを今後の研究に活 かしていきたいと考えている. 6. 謝辞 今回,筆者に武者修業の機会を与えてくださいました,21 世紀 COE プログラム「動的機能機 械システムの数理モデルと設計論」―複雑系の科学による機械工学の新たな展開― 拠点リーダー 航空宇宙工学専攻 土屋和雄教授,材料グループリーダーであり,筆者の指導教官でもある機械理 工学専攻 北村隆行教授をはじめとする関係者の皆様に厚く御礼申し上げます.また,今回の滞在 を許可していただきました,マサリック大学の Mojmir Šob 教授および研究室の皆様,ルーヴァン カトリック大学の Xavier Gonze 教授および研究室の皆様に心より感謝いたします. 参考文献 [1] X.Gonze et al., First-principles computation of material properties: the ABINIT software project, Comp.Mater.Sci., 25, (2002), 478-492. [2] X.Gonze et al., A brief introduction to the ABINIT software package, Z.Kristallogr., 220, (2005), 558-562. [3] ABINIT 公式サイト : http://www.abinit.org アムステルダムでの骨細胞武者修行報告 工学研究科機械理工学専攻 田中 基嗣 1. はじめに 私は,21 世紀 COE プログラム「動的機能機械システムの数理モデルと設計論」―複雑系の科学による 機械工学の新たな展開― における若手武者修行制度の助成を受け,2006 年 3 月 2 日から 3 月 15 日ま での約 2 週間オランダ・アムステルダムに滞在した.今回は,オランダ・アムステルダムのフライエ大 学(Vrije University)の Academic Centre for Dentistry (ACTA)の Department of Oral Cell Biology の Jenneke Klein-Nulend 教授の研究グループを訪問した.彼女の研究グループは,世界で初めて骨組織からの骨細 胞単離に成功した研究グループであり,骨細胞の活動により破骨細胞が如何に活性化されるかといった 研究テーマにおいて活発な研究活動を行っている.そのため,骨細胞研究に関して活発な討論をするこ とを目的とした.また,昨年度に彼女の研究グループを訪問した際に,ニワトリ胚から採取・単離した 骨細胞にのみ選択的に接着できる特殊な抗体(世界で彼女の研究グループだけが作ることができる!) を,私が必要になったときには譲って頂けると約束していただいた経緯がある.私の研究計画において いよいよ本格的にこの抗体が必要になってきたため,今回滞在中にその具体的な打合せをすることも重 要な目的であった.本報告では,約 2 週間の間に私が得たことについて記す. 2. フライエ大学滞在 2.1 アムステルダム(Amsterdam)とその周辺都市 オランダの首都であるアムステルダムの人口は約 75万人であり,東京駅のモデルと言われるアムステルダム中央駅(図1)を中心に市街・運河が同心円 状に広がっている.13世紀にアムステル川の河口近くにダムを造り人々が住み始めたのが「アムステル ダム」の起源であり,17世紀には海運業と貿易で栄え,黄金時代を迎えた.運河に沿って建ち並ぶ赤レ ンガや切妻屋根の建物(図2)が当時の面影を伝えており,17世紀の建物が7000以上も残っている.世 界的に有名な「アンネ・フランクの家」や「国立ゴッホ美術館」は,相変わらず世界中から来た観光客 の長蛇の列であった.オランダでは,典型的な「オランダ料理」を外食で味わうことはなかなかできな いが,いろいろな国・地域の料理を味わうことができるのが特徴である.特に,オランダの植民地であ ったインドネシアや,オランダが東インド会社を設立したことでも知られるインドの料理屋さんは非常 に多く,値段も比較的安い.昨年度のアムステルダム滞在時に「次に来たときには絶対食べに行こう!」 と心に誓ったAlbert Cuyp通り(青空日用品・食料品市で有名)近くのクルド料理屋さんに行ってみたが, すでにトルコ料理屋さんに変わっていた.時間の流れの速さを実感した次第である.しかし,すぐ近く にもう一軒別のクルド料理屋さんを見つけた. 「クルド料理」といってもどのあたりがクルド料理なのか よくわからなかったが,羊肉を使った料理がとても多いと感じた.とにかく,アムステルダムの他のレ ストランよりも安くて手が込んだ料理が出てきたことだけは確かである. 週末には,アムステルダムから国鉄に乗って,デルフト(Delft) ,デン・ハーグ(Den Haag) ,ザーンセ・ スカンス(Zaanse Schans)などに観光に行ってきた.デルフト(図3)では,その昔オランダが宗主国 図1 アムステルダム中央駅 図2 アムステルダム市街地 図3 デルフト市庁舎 図4 ビネンホフ 図7 北海 図5 国際司法裁判所 図8 風車 図6 マドゥローダム 図9 チーズ工房 スペインに対して独立戦争を行った際の指導者であり,現在のオランダ王家のご先祖様であるオレンジ 公ウィレムI世が住んでいたプリンセンホフに訪れた.1階から2階に登る階段付近に,ウィレムI世が 独立を目前にしてスペインの刺客の手に倒れた際の銃弾の痕が生々しく残っていた.デン・ハーグはオ ランダの政治の中心地であり,国会などの主要機関の入っているビネンホフと言われる歴史的な建築群 (図4) ,国連の国際司法裁判所(別名:平和宮,図5) ,マドゥローダム(世界に有名なミニチュアテ ーマパーク,図6)を訪れた.また,デン・ハーグに隣接するスエーフェニンゲン(Scheveningen)にも 訪れ,北海の冬のビーチを堪能してきた.浜辺から突き出たピア(図7)の先端には,バンジージャン プがあった.冬で非常に強い風が吹いていたのでお休みであったが,北海へのダイブはさぞかし寒いだ ろうと想像した.ザーンセ・スカンスはアムステルダムの近くにあり,日本で言うところの「明治村」 のようなものである.歴史的な民家や風車(図8)を残している地区で,チーズ工房の見学(図9)も できた. 2.2 フライエ大学 フライエ大学(図10)は,アムステルダム市の南端(ワールドトレードセンター の近く)に位置し,アムステルダム中央駅からトラム5番に乗って約 30 分程度のところにある.Jenneke の研究グループは,博士課程学生を含めて約 15 名程度のスタッフで構成されており,充分に広い居室ス ペース・実験スペースを有している.Jenneke の研究グループでは,Department of Oral Cell Biology の Head である Vincent Everts 教授の研究グループと共同で,毎週月曜日の午前 10:00 から週の始めのミーティ ングを全員で行っていた(図11) .滞在期間の間,時間を見つけては,Jenneke の研究グループを中心 に Discussion を行った.特に,Jenneke と今回も滞在のための橋渡しを行ってくれた Rommel G. Bacabac 博士(以下 Mel と呼ぶ)とは,骨細胞研究に関する今後の展望と今後共同研究を行っていくことができ るかということに関しても討論を重ねることができた.Jenneke と Mel からは, 「来年度にも京都大学に 訪問し,今後の共同研究活動に関するより具体的な討論をしたい,場合によっては早くも共同研究を開 図10 フライエ大学 図11 月曜日の定例ミーティングの様子 始したい」との期待を持っていることを伝えられた.また,共同研究をより効率よく進めるために,時 期を調整して 1 年程度留学してはどうかという話になった.1 年もあれば連名の論文を書くことも期待 でき,今後その時期と期間を調整していきたい.また,Jenneke は The City College of New York の Departments of Biomedical and Mechanical Engineering にいる Stephen C. Cowin 教授(骨細胞研究の大御所 的存在)と大変懇意にしているようで,将来的には日米欧(オランダ)を股にかけた骨細胞研究ネット ワークを形成したいと Jenneke から伝えられた. Jenneka の研究室には,Jack J.W.A. van Loon 博士というオランダでの宇宙実験・無重力状態のシミュレ ート実験等のコーディネーターが所属している.Mel も Jack とともに骨細胞株の無重力状態での応答に 関する研究を行っており,今後 Jenneke の研究グループとの共同研究が進めば,われわれの単離骨細胞 のカルシウム応答に関する実験系を無重力シミュレータに組み込んで,新しい研究テーマを共同で進め ていける可能性を感じた.また,Ton Brockers 博士,Cor Semeins 博士,Behrouz Zandieh Doulabi 博士, Marlene Knipenberg さん(博士課程学生) ,Ana Santos さん(博士課程学生) ,ZuFu Lu 氏(博士課程学生) と討論した.特に,Cor は,現在 Jenneke の研究グループで唯一骨細胞を単離し抗体を扱うことができる 人物であり,彼から抗体を使用する際のプロトコルを頂いたり,細胞に流体せん断刺激を与えるチャン バを使った実験を見学させてもらったりした.Jenneke と Cor と3人で打合わせし,私の帰国後に抗体を 使った実験計画とそのために必要な抗体の量を伝えれば,その都度抗体を空輸してもらえる算段となっ た.また,Ton,Behrouz,Marlene,Lu さんは,幹細胞を使った脊椎融合術に関する研究を行っている. ポルトガルからの留学生の Ana は博士課程を開始したばかりで,骨細胞⇔破骨細胞間の細胞コミュニケ ーションの研究に着手したばかりとのことであった. 3 月3 日の金曜日には, Ana の部屋で, Mel とMarlene と合計 4 人でメキシコ料理のデリバリーを注文してパーティをした.一軒家の最上階の部屋を間借りし ていたが,なかなかいい感じの部屋であった.アムステルダムでは Housing がもっとも難しいことだそ うで,学生には一般的な下宿の仕方のようである.3 月 7 日(15:30~16:30)には,私の研究テーマ に関するセミナを開催した.30 人程度が参加し,生体組織から単離した初代の骨細胞を使ったカルシウ ム応答(一連の細胞応答カスケードにおける最初の応答と言われている)を詳細に解明するための試み がオリジナルであると非常に好評であった. また,フライエ大学の Department of Physics の Christoph Schmidt 教授の研究グループに所属している Daisuke Mizuno 博士(Mel の共同研究者)とも有意義な討論を重ねることができ,光ピンセットとマイ クロビーズを使った細胞の刺激実験の様子を見学させてもらった.この 研究グループのボスである Christoph は現在 Nature に投稿する論文を集 中して執筆しており,議論を行う時間は残念ながらなかったが,来年に は研究グループごとドイツの University of Göttingen に移る予定とのこと で,Daisuke を介してこの研究グループとも連携することができる可能 性も感じた.他には,Vincent や Roel Breuls 博士(フライエ大学の Medical Center の Department of Physiccs and Medical Technology 所属で Mel の Supervisor の一人)とも有意義な討論を行った.いま一人,Roel の同僚 図12 Mel 料理中! であり Mel の Supervisor の一人でもある Theo H. Smit 博士と討論する時 間が残念ながら確保できなかった.次回は是非 Theo と議論をしたいと 思う. 3 月 4 日の土曜日には,Mel の部屋に遊びに行った.彼の部屋は Frederik-Hendrik 通りに位置し,日本のいわゆるワンルームマンションよ りはかなり広いワンルームのアパートである.晩ごはんにお手製のパス タをご馳走してくれた(図12) .彼の部屋にはもの凄い量の DVD があ り(図13) ,プロジェクタを使って壁に映画を映し出せるミニシアター のようになっていた. (ちなみに, パスタを食べながら少林サッカーとカ ンフーハッスルを見ました. )3 月 11 日の土曜日には,Jenneke の家(ア ンネ・フランクの家の近く)に招待されたので,Mel と二人で行った(図 図13 大量の DVD 図14 Jenneke の家で 図15 フライエ大学ワイン 14) .オランダ名物のニシンやこの日のために開封してくれたフライエ大学ワイン(図15)などをご 馳走になり,研究のこと,世界情勢のこと,日本のこと,オランダのことなど三人で夜中の二時まで話 し込んだ.もちろんバスもトラムも動いておらず,中央駅付近のホテルまで延々と歩いて帰る破目にな った.Jenneke と Mel と将来の共同研究やお互いの頻繁な行き来を堅く約束し,帰国した次第である. 3. おわりに 今回,わずか 2 週間であったが,フライエ大学に滞在することができ, 「骨細胞研究に関する討論」お よび「抗体譲渡に関する打合わせ」という目的を達成できた.また,今後の共同研究についても話を進 めることができた.今回得たことを今後の研究に活かし, 「生体組織マトリクスの変形・損傷場と細胞ネ ットワークの相互作用」に関する研究を遂行していきたいと思う. 4. 謝辞 今回,筆者に武者修行の機会を与えてくださいました,21 世紀 COE プログラム「動的機能機械シス テムの数理モデルと設計論」―複雑系の科学による機械工学の新たな展開―拠点リーダーの土屋和雄教 授(航空宇宙工学専攻) ,材料グループリーダーの北村隆行教授(機械理工学専攻)をはじめとする関係 者の皆様に厚く御礼申し上げます.また,今回の筆者の滞在を許可していただきました,フライエ大学 の Jenneke Klein-Nulend 教授,Romel G. Bacabac 氏および研究グループの皆様に心より感謝いたします.