知財求人 - 知財ポータルサイト「IP Force」
特開2022-36761シミュレーション方法、シミュレーション装置、シミュレーションプログラム、及び記録媒体
(19)【発行国】日本国特許庁(JP)
(12)【公報種別】公開特許公報(A)
(11)【公開番号】P2022036761
(43)【公開日】2022-03-08
(54)【発明の名称】シミュレーション方法、シミュレーション装置、シミュレーションプログラム、及び記録媒体
(51)【国際特許分類】
G06F 17/18 20060101AFI20220301BHJP
【ＦＩ】
G06F17/18 D
【審査請求】未請求
【請求項の数】8
【出願形態】ＯＬ
(21)【出願番号】P 2020141132
(22)【出願日】2020-08-24
【国等の委託研究の成果に係る記載事項】（出願人による申告）令和２年度、文部科学省、光・量子飛躍フラッグシッププログラム（Ｑ－ＬＥＡＰ）事業、「先端レーザーイノベーション拠点「光量子科学によるものづくりＣＰＳ化拠点」部門（光量子科学によるものづくりＣＰＳ化拠点）」委託研究、産業技術力強化法第１７条の適用を受ける特許出願
(71)【出願人】
【識別番号】301032942
【氏名又は名称】国立研究開発法人量子科学技術研究開発機構
(74)【代理人】
【識別番号】110000338
【氏名又は名称】特許業務法人ＨＡＲＡＫＥＮＺＯ ＷＯＲＬＤ ＰＡＴＥＮＴ ＆ ＴＲＡＤＥＭＡＲＫ
(72)【発明者】
【氏名】乙部 智仁
【テーマコード（参考）】
5B056
【Ｆターム（参考）】
5B056BB12
5B056BB64
(57)【要約】
【課題】変動する電磁場が印加された金属結晶中の電子系の運動を計算することが可能なシミュレーション方法を実現する。
【解決手段】シミュレーション方法（Ｓ１）は、変動する電磁場が印加された金属結晶中の電子系の運動を数値計算するシミュレーション方法である。シミュレーション方法（Ｓ１）は、少なくとも１つのプロセッサが、相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップ（Ｓ１１）を含んでいる。
【選択図】図２
【解決手段】シミュレーション方法（Ｓ１）は、変動する電磁場が印加された金属結晶中の電子系の運動を数値計算するシミュレーション方法である。シミュレーション方法（Ｓ１）は、少なくとも１つのプロセッサが、相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップ（Ｓ１１）を含んでいる。
【選択図】図２
【特許請求の範囲】
【請求項1】
変動する電磁場が印加された金属結晶中の電子系の運動を数値計算するシミュレーション方法であって、
少なくとも１つのプロセッサが、相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップを含んでいる、
ことを特徴とするシミュレーション方法。
少なくとも１つのプロセッサが、相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップを含んでいる、
ことを特徴とするシミュレーション方法。
【請求項2】
前記Vlasov方程式は、前記金属結晶中の原子・電子間の相互作用を表す結晶場ポテンシャルを含み、
当該シミュレーション方法は、前記プロセッサが、波数空間におけるAshcroft型の疑ポテンシャルの和をフーリエ変換することによって、周期境界条件を満たす前記結晶場ポテンシャルを算出するステップを更に含んでいる、
ことを特徴とする請求項１に記載のシミュレーション方法。
当該シミュレーション方法は、前記プロセッサが、波数空間におけるAshcroft型の疑ポテンシャルの和をフーリエ変換することによって、周期境界条件を満たす前記結晶場ポテンシャルを算出するステップを更に含んでいる、
ことを特徴とする請求項１に記載のシミュレーション方法。
【請求項3】
前記Vlasov方程式は、前記金属結晶中の電子・電子間の相互作用を表すハートリーポテンシャルを含み、
当該シミュレーション方法は、前記プロセッサが、波数空間における電子系の密度分布をフーリエ変換することによって、周期境界条件を満たす前記ハートリーポテンシャルを算出するステップを更に含んでいる、
ことを特徴とする請求項１又は２に記載のシミュレーション方法。
当該シミュレーション方法は、前記プロセッサが、波数空間における電子系の密度分布をフーリエ変換することによって、周期境界条件を満たす前記ハートリーポテンシャルを算出するステップを更に含んでいる、
ことを特徴とする請求項１又は２に記載のシミュレーション方法。
【請求項4】
前記プロセッサが、前記分布関数から、前記電子系の密度分布、前記電子系のエネルギー分布、前記金属結晶中の各原子の位置、前記金属結晶中の各原子のエネルギー、前記金属結晶中の電流分布、前記金属結晶中の電界分布、及び、前記金属結晶中の磁界分布の少なくとも何れかを算出するステップを更に含んでいる、
ことを特徴とする請求項１～３の何れか一項に記載のシミュレーション方法。
ことを特徴とする請求項１～３の何れか一項に記載のシミュレーション方法。
【請求項5】
前記電磁場は、レーザ光を照射したときに前記金属結晶に印加される電磁場である、
ことを特徴とする請求項１～４の何れか一項に記載のシミュレーション方法。
ことを特徴とする請求項１～４の何れか一項に記載のシミュレーション方法。
【請求項6】
変動する電磁場が印加された金属結晶中の電子系の運動を数値計算するシミュレーション装置であって、
相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップを実行する少なくとも１つのプロセッサを備えている、
ことを特徴とするシミュレーション装置。
相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップを実行する少なくとも１つのプロセッサを備えている、
ことを特徴とするシミュレーション装置。
【請求項7】
変動する電磁場が印加された金属結晶中の電子系の運動を数値計算するシミュレーションプログラムであって、
相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップを、少なくとも１つのプロセッサに実行させる、
ことを特徴とするシミュレーションプログラム。
相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップを、少なくとも１つのプロセッサに実行させる、
ことを特徴とするシミュレーションプログラム。
【請求項8】
請求項７に記載のシミュレーションプログラムが記録されているコンピュータ読み取り可能な記録媒体。
【発明の詳細な説明】
【技術分野】
【0001】
本発明は、変動する電磁場が印加された（例えば、レーザ光が照射された）金属結晶中の電子系の運動をシミュレーションする方法、装置、及びプログラムに関する。また、本発明は、そのようなプログラムが記録された記録媒体に関する。
【背景技術】
【0002】
金属中の電子系の運動を数値的に計算する方法として、時間依存密度汎関数法が知られている。しかしながら、時間依存密度汎関数法は、計算コストが高い。このため、レーザ光が照射された金属結晶中の電子系の運動など、産業上のニーズが高い現実的な電子系の運動の解析には適さない。また、より計算コストの低い方法として、時間依存密度汎関数法の基礎方式であるKohn-Sham方程式の古典極限（ｈ→０）として得られるVlasov方程式
を半古典的に解く疑似粒子解法が知られている。Vlasov方程式の疑似粒子解法を開示した文献としては、例えば、非特許文献１が挙げられる。
を半古典的に解く疑似粒子解法が知られている。Vlasov方程式の疑似粒子解法を開示した文献としては、例えば、非特許文献１が挙げられる。
【先行技術文献】
【非特許文献】
【0003】
【非特許文献1】T. Fennel, G. F. Bertsch and K. -H. Meiwes-Broer, Eur. Phys. J. D 29, 367-378 (2004)
【発明の概要】
【発明が解決しようとする課題】
【0004】
しかしながら、従来のVlasov方程式の疑似粒子解法には、孤立系にしか適用できないという制約があった。このため、例えば、レーザ光が照射された金属結晶中の電子系の運動を、従来のVlasov方程式の疑似粒子解法を用いて計算することは困難であった。
【0005】
本発明は、上記の問題に鑑みてなされたものであり、変動する電磁場が印加された金属結晶中の電子系の運動を計算することが可能なシミュレーション方法、シミュレーション装置、シミュレーションプログラムを実現することにある。
【課題を解決するための手段】
【0006】
本発明の態様１に係るシミュレーション方法は、変動する電磁場が印加された金属結晶中の電子系の運動を数値計算するシミュレーション方法であって、少なくとも１つのプロセッサが、相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップを含んでいる、という構成が採用されている。
【0007】
上記の構成によれば、疑似粒子解法を用いてVlasov方程式を周期境界条件下で解いているので、変動する電磁場が印加された金属結晶中の電子系の運動を好適に計算することができる。
【0008】
本発明の態様２に係るシミュレーション方法においては、態様１の構成に加えて、前記Vlasov方程式は、前記金属結晶中の原子・電子間の相互作用を表す結晶場ポテンシャルを含み、当該シミュレーション方法は、前記プロセッサが、波数空間におけるAshcroft型の疑ポテンシャルの和をフーリエ変換することによって、周期境界条件を満たす前記結晶場ポテンシャルを算出するステップを更に含んでいる、という構成が採用されている。
【0009】
上記の構成によれば、周期境界条件を満たす結晶場ポテンシャルを用いているので、変動する電磁場が印加された金属結晶中の電子系の運動を更に好適に計算することができる。また、特に、上記の構成によれば、Ashcroft型の疑ポテンシャルの和をフーリエ変換することによって得られた結晶場ポテンシャルを用いているので、任意の原子により構成された金属結晶中の電子系の運動を好適に計算することができる。
【0010】
本発明の態様３に係るシミュレーション方法においては、態様１又は２の構成に加えて、前記Vlasov方程式は、前記金属結晶中の電子・電子間の相互作用を表すハートリーポテンシャルを含み、当該シミュレーション方法は、前記プロセッサが、波数空間における電子系の密度分布をフーリエ変換することによって、周期境界条件を満たす前記ハートリーポテンシャルを算出するステップを更に含んでいる、という構成が採用されている。
【0011】
上記の構成によれば、周期境界条件を満たすハートリーポテンシャルを用いているので、変動する電磁場が印加された金属結晶中の電子系の運動を更に好適に計算することができる。
【0012】
本発明の態様４に係るシミュレーション方法においては、態様１～３の何れかの構成に加えて、前記プロセッサが、前記分布関数から、前記電子系の密度分布、前記電子系のエネルギー分布、前記金属結晶中の各原子の位置、前記金属結晶中の各原子のエネルギー、前記金属結晶中の電流分布、前記金属結晶中の電界分布、及び、前記金属結晶中の磁界分布の少なくとも何れかを算出するステップを更に含んでいる、という構成が採用されている。
【0013】
上記の構成によれば、前記電子系の密度分布、前記電子系のエネルギー分布、前記金属結晶中の各原子の位置、前記金属結晶中の各原子のエネルギー、前記金属結晶中の電流分布、前記金属結晶中の電界分布、及び、前記金属結晶中の磁界分布の少なくとも何れかの時間発展を好適に計算することができる。
【0014】
本発明の態様５に係るシミュレーション方法においては、態様１～４の何れかの構成に加えて、前記電磁場は、レーザ光を照射したときに前記金属結晶に印加される電磁場である、という構成が採用されている。
【0015】
上記の構成によれば、レーザが照射された金属結晶中の電子系の運動を好適に計算することができる。
【0016】
本発明の態様６に係るシミュレーション装置は、変動する電磁場が印加された金属結晶中の電子系の運動を数値計算するシミュレーション装置であって、相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップを実行する少なくとも１つのプロセッサを備えている、という構成が採用されている。
【0017】
上記の構成によれば、疑似粒子解法を用いてVlasov方程式を周期境界条件下で解いているので、変動する電磁場が印加された金属結晶中の電子系の運動を好適に計算することができる。
【0018】
本発明の態様７に係るシミュレーションプログラムは、変動する電磁場が印加された金属結晶中の電子系の運動を数値計算するシミュレーションプログラムであって、相空間における前記電子系の分布関数の時間発展を計算するために、疑似粒子解法を用いて当該分布関数が従うVlasov方程式を周期境界条件下で解くステップを、少なくとも１つのプロセッサに実行させる、という構成が採用されている。
【0019】
上記の構成によれば、前記シミュレーションプログラムを実行したコンピュータによって、変動する電磁場が印加された金属結晶中の電子系の運動を好適に計算することができる。
【0020】
本発明の態様８に係る記録媒体は、態様７のシミュレーションプログラムが記録されているコンピュータ読み取り可能な記録媒体である。
【0021】
上記の構成によれば、前記記録媒体を読み取り、前記シミュレーションプログラムを実行したコンピュータによって、変動する電磁場が印加された金属結晶中の電子系の運動を好適に計算することができる。
【発明の効果】
【0022】
本発明の一態様によれば、変動する電磁場が印加された金属結晶中の電子系の運動を好適に計算することができる。
【図面の簡単な説明】
【0023】
【発明を実施するための形態】
【0024】
【0025】
シミュレーション装置１は、シミュレーションプログラムＰ１に従って、本発明の一実施形態に係るシミュレーション方法Ｓ１を実施するための装置である。ここで、シミュレーション方法Ｓ１とは、変動する電場が印加された金属結晶（例えば、レーザ光が照射された金属結晶）中の電子系の運動を数値計算する方法のことを指す。シミュレーション方法Ｓ１の流れについては、参照する図面を代えて後述する。
【0026】
シミュレーション装置１は、汎用コンピュータを用いて実現されており、図１に示すように、プロセッサ１１と、一次メモリ１２と、二次メモリ１３と、入出力インタフェース１４と、通信インタフェース１５と、バス１６とを備えている。シミュレーション装置１において、プロセッサ１１、一次メモリ１２、二次メモリ１３、入出力インタフェース１４、及び通信インタフェース１５は、バス１６を介して相互に接続されている。
【0027】
シミュレーションプログラムＰ１は、例えば、二次メモリ１３に格納されている。この場合、シミュレーション方法Ｓ１の実施するために、プロセッサ１１は、（１）二次メモリ１３に格納されたシミュレーションプログラムＰ１を一次メモリ１２上に展開し、（２）一次メモリ１２上に展開されたシミュレーションプログラムＰ１に含まれる各命令を実
行する。なお、シミュレーションプログラムＰ１が二次メモリ１３に格納されているとは、そのソースコード、又は、そのソースコードをコンパイルすることにより得られた実行形式ファイルが二次メモリ１３に記憶されていることを指す。
行する。なお、シミュレーションプログラムＰ１が二次メモリ１３に格納されているとは、そのソースコード、又は、そのソースコードをコンパイルすることにより得られた実行形式ファイルが二次メモリ１３に記憶されていることを指す。
【0028】
プロセッサ１１として利用可能なデバイスとしては、例えば、ＣＰＵ（Central Processing Unit）、ＧＰＵ（Graphic Processing Unit）、ＤＳＰ（Digital Signal Processor）、ＭＰＵ（Micro Processing Unit）、ＦＰＵ（Floating point number Processing Unit）、ＰＰＵ（Physics Processing Unit）、マイクロコントローラ、又は、これらの組
み合わせを挙げることができる。プロセッサ１１は、「演算装置」と呼ばれることもある。
み合わせを挙げることができる。プロセッサ１１は、「演算装置」と呼ばれることもある。
【0029】
一次メモリ１２として利用可能なデバイスとしては、例えば、ＲＡＭ（Random Access Memory）を挙げることができる。一次メモリ１２は、「主記憶装置」と呼ばれることもある。二次メモリ１３として利用可能なデバイスとしては、例えば、ＨＤＤ（Hard Disk Drive）、ＳＳＤ（Solid State Drive）、ＯＤＤ（Optical Disk Drive）、フラッシュメモリ、又は、これらの組み合わせを挙げることができる。二次メモリ１３は、「補助記憶装置」と呼ばれることもある。
【0030】
入出力インタフェース１４には、入力デバイス及び／又は出力デバイスが接続される。入出力インタフェース１４としては、例えば、ＵＳＢ（Universal Serial Bus）、ＳＣＳＩ（Small Computer System Interface）、ＰＣＩ（Peripheral Component Interconnect）、HDMI（登録商標）（High Definition Multimedia Interface）、ＤＶＩ（Digital Visual Interface）などのインタフェースが挙げられる。入出力インタフェース１４に接続される入力デバイスとしては、例えば、キーボード、マウス、タッチパッド、カメラ、マイク、又は、これらの組み合わせが挙げられる。シミュレーション装置１がユーザから取得するデータは、これらの入力デバイスを介してシミュレーション装置１に入力され、例えば、一次メモリ１２に記憶される。入出力インタフェース１４に接続される出力デバイスとしては、例えば、ディスプレイ、プロジェクタ、スピーカ、ヘッドホン、プリンタ、又は、これらの組み合わせが挙げられる。シミュレーション装置１がユーザに提供する情報は、これらの出力デバイスを介してシミュレーション装置１から出力される。
【0031】
なお、シミュレーション装置１は、ラップトップ型コンピュータのように、入力デバイスとして機能するキーボードと、出力デバイスとして機能するディスプレイとを、それぞれ内蔵してもよい。或いは、シミュレーション装置１は、タブレット型コンピュータのように、入力デバイス及び出力デバイスの両方として機能するタッチパネルを内蔵していてもよい。
【0032】
通信インタフェース１５には、通信ネットワークを介して他のコンピュータが有線接続又は無線接続される。通信インタフェース１５としては、例えば、イーサネット（登録商標）、Ｗｉ－Ｆｉ（登録商標）などが挙げられる。利用可能な通信ネットワークとしては、例えば、ＰＡＮ（Personal Area Network）、ＬＡＮ（Local Area Network）、ＣＡＮ
（Campus Area Network）、ＭＡＮ（Metropolitan Area Network）、ＷＡＮ（Wide Area Network）、ＧＡＮ（Global Area Network）、又は、これらの通信ネットワークを含むインターネットワークが挙げられる。インターネットワークは、イントラネットであってもよいし、エクストラネットであってもよいし、インターネットであってもよい。シミュレーション方法Ｓ１においてシミュレーション装置１が他のコンピュータから取得するデータ、及び、シミュレーション方法Ｓ１においてシミュレーション装置１が他のコンピュータに提供するデータは、これらの通信ネットワークを介して送受信される。
（Campus Area Network）、ＭＡＮ（Metropolitan Area Network）、ＷＡＮ（Wide Area Network）、ＧＡＮ（Global Area Network）、又は、これらの通信ネットワークを含むインターネットワークが挙げられる。インターネットワークは、イントラネットであってもよいし、エクストラネットであってもよいし、インターネットであってもよい。シミュレーション方法Ｓ１においてシミュレーション装置１が他のコンピュータから取得するデータ、及び、シミュレーション方法Ｓ１においてシミュレーション装置１が他のコンピュータに提供するデータは、これらの通信ネットワークを介して送受信される。
【0033】
なお、本実施形態においては、単一のプロセッサ（プロセッサ１１）を用いてシミュレ
ーション方法Ｓ１を実施する構成を採用しているが、本発明は、これに限定されない。すなわち、複数のプロセッサを用いてシミュレーション方法Ｓ１を実施する構成を採用してもよい。この場合、シミュレーション方法Ｓ１を実施する複数のプロセッサは、単一のコンピュータに集中して設けられ、バスを介して相互に通信可能に構成されていてもよい。或いは、シミュレーション方法Ｓ１を実施する複数のプロセッサは、複数のコンピュータに分散して設けられ、ネットワークを介して相互に通信可能に構成されていてもよい。
ーション方法Ｓ１を実施する構成を採用しているが、本発明は、これに限定されない。すなわち、複数のプロセッサを用いてシミュレーション方法Ｓ１を実施する構成を採用してもよい。この場合、シミュレーション方法Ｓ１を実施する複数のプロセッサは、単一のコンピュータに集中して設けられ、バスを介して相互に通信可能に構成されていてもよい。或いは、シミュレーション方法Ｓ１を実施する複数のプロセッサは、複数のコンピュータに分散して設けられ、ネットワークを介して相互に通信可能に構成されていてもよい。
【0034】
また、本実施形態においては、二次メモリ１３がシミュレーション装置１に内蔵されている構成が採用されているが、本発明は、これに限定されない。すなわち、二次メモリ１３が通信インタフェース１５を介してシミュレーション装置１と接続された他のコンピュータに内蔵されている構成を採用してもよい。或いは、二次メモリ１３が入出力インタフェース１４を介してシミュレーション装置１と接続されたドライブ装置に内蔵されている構成を採用してもよい。
【0035】
また、本実施形態においては、シミュレーション装置１における記憶を２つのメモリ（一次メモリ１２及び二次メモリ１３）により実現する構成を採用しているが、本発明は、これに限定されない。すなわち、シミュレーション装置１における記憶を１つのメモリにより実現する構成を採用してもよい。この場合、例えば、そのメモリの或る記憶領域を一次メモリ１２として利用し、そのメモリの他の記憶領域を二次メモリ１３として利用すればよい。
【0036】
なお、シミュレーションプログラムＰ１は、任意の記録媒体、特に、コンピュータ読み取り可能な、有形の、一時的でない記録媒体に記録され得る。コンピュータ読み取り可能な、有形の、一時的でない記録媒体としては、上述した二次メモリ１３の他に、例えば、磁気テープ、磁気ディスク、光ディスク、光磁気ディスク、フラッシュメモリなどが挙げられる。また、シミュレーションプログラムＰ１は、任意の伝送媒体を介して伝送され得る。伝送媒体としては、上述した通信ネットワークの他に、例えば放送波などが挙げられる。
【0037】
（Vlasov方程式の疑似粒子解法）
変動する電磁場が印加された金属結晶中の電子の運動は、Kohn-Sham方程式により記述
される。Vlasov方程式は、Kohn-Sham方程式の古典極限（ｈ→０）として得られる、相空
間における電子系の分布関数ｆ（ｒ，ｐ，ｔ）が従う方程式である。Vlasov方程式は、以下のように書き下すことができる。
変動する電磁場が印加された金属結晶中の電子の運動は、Kohn-Sham方程式により記述
される。Vlasov方程式は、Kohn-Sham方程式の古典極限（ｈ→０）として得られる、相空
間における電子系の分布関数ｆ（ｒ，ｐ，ｔ）が従う方程式である。Vlasov方程式は、以下のように書き下すことができる。
【数1】
【0038】
ここで、ｒは位置、ｐは運動量、ｔは時間、ｍは電子の質量である。また、ＩUUは、衝突項と呼ばれる定数である。電子間衝突を考慮する場合には、ＩUUとして分布関数ｆ（ｒ，ｐ，ｔ）の汎関数ＩUU＝Ｄt［ｆ（ｒ，ｐ，ｔ）］を用い（Physical Review Letters 81, 5524(1998)参照）、電子間衝突を無視する場合には、ＩUU＝０とする。また、Ｖps（
ｒ）は、電子－原子（イオン）間の相互作用を表す結晶場ポテンシャルであり、ＶH［ρ
（ｒ，ｔ）］は、電子－電子間の相互作用を表すハートリーポテンシャルであり、Ｖxc［
ρ（ｒ，ｔ）］は、局所電子密度による相関交換ポテンシャルであり、Ｖext（ｒ，ｔ）
は、金属結晶に印加される電磁場による外部ポテンシャルである。ここで、ρ（ｒ，ｔ）は、電子系の密度分布であり、ρ（ｒ，ｔ）＝∫ｄｐｆ（ｒ，ｐ，ｔ）により定義される。結晶場ポテンシャルＶps（ｒ）とハートリーポテンシャルＶH［ρ（ｒ，ｔ）］と相関
交換ポテンシャルＶxc［ρ（ｒ，ｔ）］との和を、実効ポテンシャルＶeff（ｒ，ｔ）と
呼び、実効ポテンシャルＶeff（ｒ，ｔ）と外部ポテンシャルＶext（ｒ，ｔ）との和を、トータルポテンシャルＶtot（ｒ，ｔ）と呼ぶ。
ｒ）は、電子－原子（イオン）間の相互作用を表す結晶場ポテンシャルであり、ＶH［ρ
（ｒ，ｔ）］は、電子－電子間の相互作用を表すハートリーポテンシャルであり、Ｖxc［
ρ（ｒ，ｔ）］は、局所電子密度による相関交換ポテンシャルであり、Ｖext（ｒ，ｔ）
は、金属結晶に印加される電磁場による外部ポテンシャルである。ここで、ρ（ｒ，ｔ）は、電子系の密度分布であり、ρ（ｒ，ｔ）＝∫ｄｐｆ（ｒ，ｐ，ｔ）により定義される。結晶場ポテンシャルＶps（ｒ）とハートリーポテンシャルＶH［ρ（ｒ，ｔ）］と相関
交換ポテンシャルＶxc［ρ（ｒ，ｔ）］との和を、実効ポテンシャルＶeff（ｒ，ｔ）と
呼び、実効ポテンシャルＶeff（ｒ，ｔ）と外部ポテンシャルＶext（ｒ，ｔ）との和を、トータルポテンシャルＶtot（ｒ，ｔ）と呼ぶ。
【0039】
相空間における電子系の分布関数ｆ（ｒ，ｐ，ｔ）は、疑似粒子の集合とみなすことによって、下記の式（２）のように近似することができる。
【数2】
【0040】
ここで、Ｎppは疑似粒子の個数であり、Ｎsは疑似粒子の個数Ｎppを電子の個数Ｎeで除した商Ｎpp／Ｎeである。また、πは円周率であり、σr及びσpは幅パラメータ（Width Parameter）と呼ばれる定数である。また、ｒiはｉ番目の疑似粒子の座標であり、ｐiはｉ番目の疑似粒子の運動量である。
【0041】
Vlasov方程式（１）を上記の式（２）で近似すると共に、外部ポテンシャルＶext（ｒ
，ｔ）による力を局所的一様電場から受ける力ｅＦ（ｒ，ｔ）で近似すると、疑似粒子が従う下記のNewton方程式（３）が得られる。
，ｔ）による力を局所的一様電場から受ける力ｅＦ（ｒ，ｔ）で近似すると、疑似粒子が従う下記のNewton方程式（３）が得られる。
【数3】
【0042】
上記の式（２）を用いることによって、時刻ｔでの相空間における電子系の分布関数ｆ（ｒ，ｐ，ｔ）から、各疑似粒子の位置ｒi（ｔ）及び運動量ｐi（ｔ）を求めることができる。そして、Newton方程式（３）を解くことによって、時刻ｔでの各疑似粒子の位置ｒi（ｔ）及び運動量ｐi（ｔ）から、時刻ｔ＋Δｔでの各疑似粒子の位置ｒi（ｔ＋Δｔ）
及び運動量ｐi（ｔ＋Δｔ）を求めることができる。そして、上記の式（２）を用いるこ
とによって、時刻ｔ＋Δｔでの各疑似粒子の位置ｒi（ｔ＋Δｔ）及び運動量ｐi（ｔ＋Δｔ）から、時刻ｔ＋Δｔでの相空間における電子系の分布関数ｆ（ｒ，ｐ，ｔ＋Δｔ）を求めることができる。このように、Newton方程式（３）を解くことによって、相空間における電子系の分布関数ｆ（ｒ，ｐ，ｔ）の時間発展を求めることを、Vlasov方程式（１）の疑似粒子解法又は半古典解法という。なお、以下の説明においては、上記の式（２）のことを、変換式（２）と呼ぶことにする。
及び運動量ｐi（ｔ＋Δｔ）を求めることができる。そして、上記の式（２）を用いるこ
とによって、時刻ｔ＋Δｔでの各疑似粒子の位置ｒi（ｔ＋Δｔ）及び運動量ｐi（ｔ＋Δｔ）から、時刻ｔ＋Δｔでの相空間における電子系の分布関数ｆ（ｒ，ｐ，ｔ＋Δｔ）を求めることができる。このように、Newton方程式（３）を解くことによって、相空間における電子系の分布関数ｆ（ｒ，ｐ，ｔ）の時間発展を求めることを、Vlasov方程式（１）の疑似粒子解法又は半古典解法という。なお、以下の説明においては、上記の式（２）のことを、変換式（２）と呼ぶことにする。
【0043】
【0044】
シミュレーション方法Ｓ１は、変動する電場が印加された金属結晶（例えば、レーザ光
が照射された金属結晶）中の電子系の運動を、上述したVlasov方程式（１）の疑似粒子解法を用いて数値計算する方法である。シミュレーション方法Ｓ１は、図２に示すように、定数設定ステップＳ１１、結晶場ポテンシャル設定ステップＳ１２、初期状態設定ステップＳ１３、ハートリーポテンシャル設定ステップＳ１４、相関交換ポテンシャル設定ステップＳ１５、時間発展計算ステップＳ１６、及び物理量算出ステップＳ１７を含んでいる。ハートリーポテンシャル設定ステップＳ１４、相関交換ポテンシャル設定ステップＳ１５、時間発展計算ステップＳ１６、及び物理量算出ステップＳ１７は、繰り返し実行される。
が照射された金属結晶）中の電子系の運動を、上述したVlasov方程式（１）の疑似粒子解法を用いて数値計算する方法である。シミュレーション方法Ｓ１は、図２に示すように、定数設定ステップＳ１１、結晶場ポテンシャル設定ステップＳ１２、初期状態設定ステップＳ１３、ハートリーポテンシャル設定ステップＳ１４、相関交換ポテンシャル設定ステップＳ１５、時間発展計算ステップＳ１６、及び物理量算出ステップＳ１７を含んでいる。ハートリーポテンシャル設定ステップＳ１４、相関交換ポテンシャル設定ステップＳ１５、時間発展計算ステップＳ１６、及び物理量算出ステップＳ１７は、繰り返し実行される。
【0045】
定数設定ステップＳ１１は、Vlasov方程式（１）の疑似粒子解法において用いられる各種定数の値を、プロセッサ１１が設定するステップである。定数設定ステップＳ１１において値が設定される定数としては、前述した変換式（２）及びNewton方程式（３）に含まれる定数の他に、後述する疑ポテンシャルｗ（Ｑ）の定義式（４）、結晶場ポテンシャルＶps（ｒ）の定義式（５）、及びハートリーポテンシャルＶH［ρ（ｒ，ｔ）］の定義式
（７）に含まれる定数が挙げられる。これらの定数の値は、予め定められた値であってもよいし、入力デバイスを用いてユーザが入力した値であってもよい。
（７）に含まれる定数が挙げられる。これらの定数の値は、予め定められた値であってもよいし、入力デバイスを用いてユーザが入力した値であってもよい。
【0046】
結晶場ポテンシャル設定ステップＳ１２は、Vlasov方程式（１）に現れる結晶場ポテンシャルＶps（ｒ）を、プロセッサ１１が設定するステップである。本実施形態において、プロセッサ１１は、波数空間ＱにおけるAshcroft型の疑ポテンシャルｗ（Ｑ）の和をフーリエ変換することによって、周期境界条件Ｖps（ｒ＋Ｌ）＝Ｖps（ｒ）を満たす結晶場ポテンシャルＶps（ｒ）を算出する。一例として、プロセッサ１１は、下記の式（４）により定義される疑ポテンシャルｗ（Ｑ）の和をフーリエ変換することによって、下記の式（５）により定義される結晶場ポテンシャルＶps（ｒ）を算出する。なお、フーリエ変換を数値的に行う方法は公知であるため、ここではその詳細な説明を割愛する。
【数4】
【数5】
【0047】
上記の式（４）において、α、β＝（α３－２α）／［４（α２－１）］、Ａ＝α２／２－αβ、Ｒは、金属結晶を構成する原子に応じて定まる定数である。例えば、金属結晶を構成する原子がＡｌである場合、α＝３．６３５又は３．５７３、Ｒ＝０．３３４又は０．３１７であり、金属結晶を構成する原子がＭｇである場合、α＝３．５０２又は３．５０８であり、Ｒ＝０．３８２又は０．３８３である。また、上記の式（４）において、ｚは金属結晶を構成する原子の原子価である。また、上記の式（５）において、Ｉiはｉ
番目の原子の座標を表す定数である。また、上記の式（５）において、Ｖcellは単位胞の体積を表す定数である。
番目の原子の座標を表す定数である。また、上記の式（５）において、Ｖcellは単位胞の体積を表す定数である。
【0048】
初期状態設定ステップＳ１３は、ｔ＝０での相空間における電子系の分布関数ｆ（ｒ，ｐ，０）、ｔ＝０での疑似粒子の位置ｒi（０）及び運動量ｐi（０）、及び、ｔ＝０での電子系の密度分布ρ（ｒ，０）を、プロセッサ１１が設定するステップである。本実施形態において、プロセッサ１１は、ｔ＝０での相空間における電子系の分布関数ｆ（ｒ，ｐ，０）として、トーマス・フェルミモデルを用いて底状態に対応する電子系の分布関数を
算出する。なお、トーマス・フェルミモデルを用いて基底状態に対応する電子系の分布関数を算出する方法は公知であるため、ここではその詳細な説明を割愛する。また、本実施形態において、プロセッサ１１は、変換式（２）を用いて、ｔ＝０での相空間における電子系の分布関数ｆ（ｒ，ｐ，０）から、ｔ＝０での疑似粒子の位置ｒi（０）及び運動量
ｐi（０）を算出する。また、本実施形態において、プロセッサ１１は、ρ（ｒ，０）＝
∫ｄｐｆ（ｒ，ｐ，０）を用いて、ｔ＝０での相空間における電子系の分布関数ｆ（ｒ，ｐ，０）から、ｔ＝０での電子系の密度分布ρ（ｒ，０）を算出する。
算出する。なお、トーマス・フェルミモデルを用いて基底状態に対応する電子系の分布関数を算出する方法は公知であるため、ここではその詳細な説明を割愛する。また、本実施形態において、プロセッサ１１は、変換式（２）を用いて、ｔ＝０での相空間における電子系の分布関数ｆ（ｒ，ｐ，０）から、ｔ＝０での疑似粒子の位置ｒi（０）及び運動量
ｐi（０）を算出する。また、本実施形態において、プロセッサ１１は、ρ（ｒ，０）＝
∫ｄｐｆ（ｒ，ｐ，０）を用いて、ｔ＝０での相空間における電子系の分布関数ｆ（ｒ，ｐ，０）から、ｔ＝０での電子系の密度分布ρ（ｒ，０）を算出する。
【0049】
ハートリーポテンシャル設定ステップＳ１４は、Vlasov方程式（１）に現れるハートリーポテンシャルＶH［ρ（ｒ，ｔ）］を、プロセッサ１１が設定するステップである。本
実施形態において、プロセッサ１１は、初期状態設定ステップＳ１３又は物理量算出ステップＳ１７にて算出された配位空間における電子系の密度分布ρ（ｒ，ｔ）を逆フーリエ変換することによって、波数空間における電子系の密度分布ρG（Ｑ，ｔ）を算出する。
一例として、プロセッサ１１は、配位空間における電子系の密度分布ρ（ｒ，ｔ）を逆フーリエ変換することによって、下記の式（６）により定義される波数空間における電子系の密度分布ρG（Ｑ，ｔ）を算出する。また、本実施形態において、プロセッサ１１は、
波数空間Ｑにおける電子系の密度分布ρG（Ｑ，ｔ）をフーリエ変換することによって、
周期境界条件ＶH［ρ（ｒ＋Ｌ，ｔ）］＝ＶH［ρ（ｒ，ｔ）］を満たすハートリーポテンシャルＶH［ρ（ｒ，ｔ）］を算出する。一例として、プロセッサ１１は、波数空間にお
ける電子系の密度分布ρG（Ｑ，ｔ）をフーリエ変換することによって、下記の式（７）
により定義されるハートリーポテンシャルＶH［ρ（ｒ，ｔ）］を算出する。なお、逆フ
ーリエ変換及びフーリエ変換を数値的に行う方法は公知であるため、ここではその詳細な説明を割愛する。
実施形態において、プロセッサ１１は、初期状態設定ステップＳ１３又は物理量算出ステップＳ１７にて算出された配位空間における電子系の密度分布ρ（ｒ，ｔ）を逆フーリエ変換することによって、波数空間における電子系の密度分布ρG（Ｑ，ｔ）を算出する。
一例として、プロセッサ１１は、配位空間における電子系の密度分布ρ（ｒ，ｔ）を逆フーリエ変換することによって、下記の式（６）により定義される波数空間における電子系の密度分布ρG（Ｑ，ｔ）を算出する。また、本実施形態において、プロセッサ１１は、
波数空間Ｑにおける電子系の密度分布ρG（Ｑ，ｔ）をフーリエ変換することによって、
周期境界条件ＶH［ρ（ｒ＋Ｌ，ｔ）］＝ＶH［ρ（ｒ，ｔ）］を満たすハートリーポテンシャルＶH［ρ（ｒ，ｔ）］を算出する。一例として、プロセッサ１１は、波数空間にお
ける電子系の密度分布ρG（Ｑ，ｔ）をフーリエ変換することによって、下記の式（７）
により定義されるハートリーポテンシャルＶH［ρ（ｒ，ｔ）］を算出する。なお、逆フ
ーリエ変換及びフーリエ変換を数値的に行う方法は公知であるため、ここではその詳細な説明を割愛する。
【数6】
【数7】
【0050】
相関交換ポテンシャル設定ステップＳ１５は、Vlasov方程式（１）に現れる相関交換ポテンシャルＶxc［ρ（ｒ，ｔ）］を、プロセッサ１１が設定するステップである。本実施形態において、プロセッサ１１は、ＬＤＡ近似を用いて、初期状態設定ステップＳ１３又は物理量算出ステップＳ１７にて算出された電子密度分布ρ（ｒ，ｔ）から、相関交換ポテンシャルＶxc［ρ（ｒ，ｔ）］を算出する。ＬＤＡ近似を用いて電子密度分布ρ（ｒ，ｔ）から相関交換ポテンシャルＶxc［ρ（ｒ，ｔ）］を算出する方法は公知であるため、ここではその詳細な説明を割愛する。
【0051】
時間発展計算ステップＳ１６は、プロセッサ１１が、疑似粒子解法を用いて周期境界条件下でVlasov方程式（１）を解くステップである。本実施形態において、プロセッサ１１は、Newton方程式（３）を用いて、時刻ｔでの各疑似粒子の位置ｒi（ｔ）及び運動量ｐi（ｔ）から、時刻ｔ＋Δｔでの各疑似粒子の位置ｒi（ｔ＋Δｔ）及び運動量ｐi（ｔ＋Δｔ）を算出する。また、本実施形態において、プロセッサ１１は、変換式（２）を用いて、時刻ｔ＋Δｔでの各疑似粒子の位置ｒi（ｔ＋Δｔ）及び運動量ｐi（ｔ＋Δｔ）から、時刻ｔ＋Δｔでの相空間における電子系の分布関数ｆ（ｒ，ｐ，ｔ＋Δｔ）を算出する。なお、周期境界条件としてｆ（ｒ＋Ｌ，ｐ、ｔ）＝ｆ（ｒ，ｐ，ｔ）を課す場合、Newton方程式（３）を解く際に、位置ｒi＋ｎ×Ｌ（ｎは任意の整数）は位置ｒiと同一視される
。なお、Newton方程式（３）を用いて、時刻ｔでの各疑似粒子の位置ｒi（ｔ）及び運動
量ｐi（ｔ）から、時刻ｔ＋Δｔでの各疑似粒子の位置ｒi（ｔ＋Δｔ）及び運動量ｐi（
ｔ＋Δｔ）を数値的に算出する方法は公知であるため、ここではその詳細な説明を割愛する。
。なお、Newton方程式（３）を用いて、時刻ｔでの各疑似粒子の位置ｒi（ｔ）及び運動
量ｐi（ｔ）から、時刻ｔ＋Δｔでの各疑似粒子の位置ｒi（ｔ＋Δｔ）及び運動量ｐi（
ｔ＋Δｔ）を数値的に算出する方法は公知であるため、ここではその詳細な説明を割愛する。
【0052】
なお、Newton方程式（３）を数値的に解く方法は、上記の方法に限定されない。例えば、Newton方程式（３）は、以下のように解いてもよい。まず、時刻ｔでの各疑似粒子の位置ｒi（ｔ）及び運動量ｐi（ｔ）から、Newton方程式（３）を用いて、時刻ｔ＋Δｔ／２での各疑似粒子の位置ｒi（ｔ＋Δｔ／２）する。次に、時刻ｔ＋Δｔ／２での各疑似粒
子の位置ｒi（ｔ＋Δｔ／２）から、上述した方法を用いて、ハートリーポテンシャルＶH［ρ（ｒ，ｔ＋Δｔ）］及び相関交換ポテンシャルＶxc［ρ（ｒ，ｔ＋Δｔ）］を算出する。そして、これらポテンシャルを含むNewton方程式（３）を用いて、時刻ｔ＋Δｔでの各疑似粒子の位置ｒi（ｔ＋Δｔ）及び運動量ｐi（ｔ＋Δｔ）を算出する。
子の位置ｒi（ｔ＋Δｔ／２）から、上述した方法を用いて、ハートリーポテンシャルＶH［ρ（ｒ，ｔ＋Δｔ）］及び相関交換ポテンシャルＶxc［ρ（ｒ，ｔ＋Δｔ）］を算出する。そして、これらポテンシャルを含むNewton方程式（３）を用いて、時刻ｔ＋Δｔでの各疑似粒子の位置ｒi（ｔ＋Δｔ）及び運動量ｐi（ｔ＋Δｔ）を算出する。
【0053】
物理量算出ステップＳ１７は、時刻ｔ＋Δｔにおける各種物理量を算出するステップである。本実施形態において、プロセッサ１１は、時間発展計算ステップＳ１６にて得られた時刻ｔ＋Δｔでの電子系の分布関数ｆ（ｒ，ｐ，ｔ＋Δｔ）から、時刻ｔ＋Δｔでの電子系の密度分布ρ（ｒ，ｔ＋Δｔ）その他の物理量を算出する。物理量算出ステップＳ１７にて算出可能な、電子系の密度分布ρ（ｒ，ｔ＋Δｔ）以外の物理量としては、電子系のエネルギー分布、金属結晶を構成する各原子の位置、金属結晶を構成する各原子のエネルギー、金属結晶中の電流分布、金属結晶中の電界分布（電場）、金属結晶中の磁界分布（磁場）などが挙げられる。
【0054】
（実施例１）
レーザ光が照射されたアルミニウム単結晶中の電子系の運動を上述したシミュレーション方法Ｓ１（実施例）及び時間依存密度汎関数法（比較例）を用いて計算し、電流密度の時間変化を算出した。レーザ光については、ピーク強度を５×１０１２Ｗ／ｃｍ２、波長を４００ｎｍ、パルス半値全幅を５ｆｓに設定した。
レーザ光が照射されたアルミニウム単結晶中の電子系の運動を上述したシミュレーション方法Ｓ１（実施例）及び時間依存密度汎関数法（比較例）を用いて計算し、電流密度の時間変化を算出した。レーザ光については、ピーク強度を５×１０１２Ｗ／ｃｍ２、波長を４００ｎｍ、パルス半値全幅を５ｆｓに設定した。
【0055】
その計算結果を図３に示す。図３において、点線は上述したシミュレーション方法Ｓ１を用いて算出した電流密度であり、実線は時間依存密度汎関数法を用いて算出した電流密度である（実線のみに見える部分は点線と実線とが重なっている点に留意されたい）。上述したシミュレーション方法Ｓ１を用いて算出した電流密度の時間変化は、時間依存密度汎関数法を用いて算出した電流密度の時間変化と良く一致している。このことは、上述したシミュレーション方法Ｓ１を用いれば、時間依存密度汎関数法に匹敵する精度のシミュレーションが可能であることを示唆している。
【0056】
（実施例２）
上述したシミュレーション方法Ｓ１は、周期境界条件の周期Ｌを金属薄膜の厚みよりも大きく設定することにより、金属薄膜中の電子系の運動の計算にも利用することができる。そこで、レーザ光が照射されたアルミニウム薄膜中の電子系の運動を上述したシミュレーション方法Ｓ１を用いて計算し、電界分布、エネルギー密度、及び電子密度を算出した。アルミニウム薄膜については、厚さを０．４μｍに設定した。レーザ光については、ピーク強度を１×１０１２Ｗ／ｃｍ２、波長を４００ｎｍ、パルス半値全幅を７ｆｓに設定した。
上述したシミュレーション方法Ｓ１は、周期境界条件の周期Ｌを金属薄膜の厚みよりも大きく設定することにより、金属薄膜中の電子系の運動の計算にも利用することができる。そこで、レーザ光が照射されたアルミニウム薄膜中の電子系の運動を上述したシミュレーション方法Ｓ１を用いて計算し、電界分布、エネルギー密度、及び電子密度を算出した。アルミニウム薄膜については、厚さを０．４μｍに設定した。レーザ光については、ピーク強度を１×１０１２Ｗ／ｃｍ２、波長を４００ｎｍ、パルス半値全幅を７ｆｓに設定した。
【0057】
その計算結果を図４に示す。図４において、（ａ）は電界分布を表すグラフであり、（ｂ）は電子系のエネルギー分布を表すグラフであり、（ｃ）は電子系の密度分布を表すグラフである。各グラフにおいて、実線はパルス通過前（ｔ＝０ｆｓ）の分布、点線はパルス通過中（ｔ＝４．８４ｆｓ）の分布、鎖線はパルス通過後（ｔ＝１２．１０ｆｓ）の分
布を表す。
布を表す。
【0058】
図４の（ａ）によれば、アルミニウム薄膜に照射されたレーザ光が、アルミニウム薄膜において反射波と透過波とに分かれる様子が見て取れる。また、図４の（ｂ）によれば、アルミニウム薄膜の表面（レーザ光が入射する方の面）に加えて、アルミニウム薄膜の裏面（レーザ光が出射する方の面）及び中心部においても、電子励起が起こっていることが見て取れる。また、図４の（ｂ）によれば、アルミニウム薄膜の表面及び裏面から電子が放出されている様子が見て取れる。上述したシミュレーション方法Ｓ１によれば、このように、レーザ光が照射された金属薄膜において生じ得る様々な物理現象を、理解又は発見することができる。
【0059】
〔付記事項〕
本発明は、上述した実施形態に限定されるものでなく、請求項に示した範囲で種々の変更が可能であり、上述した実施形態にそれぞれ開示された技術的手段を適宜組み合わせて得られる実施形態についても、本発明の技術的範囲に含まれる。
本発明は、上述した実施形態に限定されるものでなく、請求項に示した範囲で種々の変更が可能であり、上述した実施形態にそれぞれ開示された技術的手段を適宜組み合わせて得られる実施形態についても、本発明の技術的範囲に含まれる。
【符号の説明】
【0060】
１ シミュレーション装置
１１ プロセッサ
１２ 一次メモリ
１３ 二次メモリ
１４ 入出力インタフェース
１５ 通信インタフェース
１６ バス
Ｓ１ シミュレーション方法
Ｓ１１ 定数設定ステップ
Ｓ１２ 結晶場ポテンシャル設定ステップ
Ｓ１３ 初期状態設定ステップ
Ｓ１４ ハートリーポテンシャル設定ステップ
Ｓ１５ 相関交換ポテンシャル設定ステップ
Ｓ１６ 時間発展計算ステップ
Ｓ１７ 物理量算出ステップ
１１ プロセッサ
１２ 一次メモリ
１３ 二次メモリ
１４ 入出力インタフェース
１５ 通信インタフェース
１６ バス
Ｓ１ シミュレーション方法
Ｓ１１ 定数設定ステップ
Ｓ１２ 結晶場ポテンシャル設定ステップ
Ｓ１３ 初期状態設定ステップ
Ｓ１４ ハートリーポテンシャル設定ステップ
Ｓ１５ 相関交換ポテンシャル設定ステップ
Ｓ１６ 時間発展計算ステップ
Ｓ１７ 物理量算出ステップ
【図1】
【図2】
【図3】
【図4】