JP2022021390A - 分析システム、情報処理装置、分析方法及びプログラム - Google Patents

Abstract

【課題】試料の結晶内部の３次元的な性質を効率的に可視化することが可能な分析システムを提供する。【解決手段】分析装置制御部１１０は、分析装置１０を制御して、位置ｋに、フェムト秒レーザパルス２を照射させる。分析装置制御部１１０は、光学遅延ステージ４０を制御して、トリガパルス６が光伝導スイッチ３２に到達する時間を制御する。これにより、位置ｋで発生したテラヘルツ波８０の時間波形が得られる。情報処理装置１００の取得部１２０は、観測されたテラヘルツ波８０を示す観測データを取得する。分析部１３０は、得られた観測データに対して上述した因子分析を行う。【選択図】図４

Description

本発明は、分析システム、情報処理装置、分析方法及びプログラムに関する。

試料の結晶内部の３次元的な性質を可視化する技術（ドメインイメージング）が提案されている。例えば、非特許文献１～４は、第二高調波発生（second-harmonic generation：ＳＨＧ）による３次元ドメイン可視化手法を開示している。また、非特許文献５は、フェムト秒レーザを結晶に照射することにより結晶から放射されるテラヘルツ（ＴＨｚ）波を観測することによってドメインを可視化する方法を開示している。また、非特許文献６は、歯にテラヘルツ波を当てて、エコーの位置からその断面構造を解析する研究について開示している。

Kaneshiro, Junichi, et al.、"Three-dimensional observations of polar domain structures using a confocal second-harmonic generation interference microscope"、Journal of Applied Physics 104.5 (2008): 054112 Sheng, Y., Best, A., Butt, H. J., Krolikowski, W., Arie, A., & Koynov, K.、"Three-dimensional ferroelectric domain visualization by Cerenkov-type second harmonic generation"、Optics express, 18(16), 16539-16545 (2010) Potnis, P. R., Tsou, N. T., & Huber, J. E.、"A review of domain modelling and domain imaging techniques in ferroelectric crystals"、Materials, 4(2), 417-447 (2011) Ayoub, M., Futterlieb, H., Imbrock, J., & Denz, C.、"3D imaging of ferroelectric kinetics during electrically driven switching"、Advanced Materials, 29(5), 1603325 (2017) M. Sotome, N. Kida, S. Horiuchi and H. Okamoto、"Terahertz radiation imaging of ferroelectric domain topography in room-temperature hydrogen-bonded supramolecular ferroelectrics"、ACS Photonics 2, 1373: 1-11 (2015) David A. Crawley, Christopher Longbottom, Vincent P. Wallace, Bryan E. Cole, Donald D. Arnone, and Michael Pepper、"Three-dimensional terahertz pulse imaging of dentaltissue"、Journal of Biomedical Optics 8(2), 303-307 (April 2003)

上述した非特許文献にかかる技術では、ドメインイメージングを行う際に、対象物質について、屈折率スペクトル等の光学スペクトル情報（分光データ）を予め準備しておく必要がある可能性がある。一方、未知の対象物質の場合、光学スペクトル情報はない。また、物質の表面粗さや形状や加工の制限のために予め光学スペクトルの取得することは困難であるおそれがある。したがって、非特許文献にかかる技術では、試料の結晶内部の３次元的な性質（ドメイン）を効率的に可視化することができないおそれがある。

本発明は、このような課題を解決するためになされたものであり、試料の結晶内部の３次元的な性質を効率的に可視化することが可能な分析システム、情報処理装置、分析方法及びプログラムを提供することを目的とする。

本発明にかかる分析システムは、試料の第１の面に沿ってレーザ光を走査して前記レーザ光を前記第１の面の側から前記試料に入射させ、前記試料から出射したテラヘルツ波を検出する分析装置と、前記テラヘルツ波の観測された振幅の時間的な推移を示す観測データを取得する取得部と、取得された前記観測データに対して統計的データ解析を行うことによって、前記試料の、前記レーザ光の入射方向に沿った方向である深さ方向の複数の深さ位置それぞれにおける、前記レーザ光を走査した位置に対応する照射位置に関するドメインの分布を推定する分析部と、を有する。

好ましくは、前記取得部は、前記試料の前記複数の深さ位置それぞれで発生し前記試料内を伝搬した後に前記試料の外に出射された前記テラヘルツ波に関する前記観測データを取得し、前記分析部は、取得された前記観測データに対して統計的データ解析を行うことによって、前記取得部で取得された前記観測データに対応する前記テラヘルツ波の波形を推定する。

好ましくは、前記分析部は、前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する位相スペクトルを用いて、推定された前記ドメインの分布に対応する前記深さ位置を推定する。

好ましくは、前記分析部は、取得された前記観測データに対して因子分析を行うことによって、前記ドメインの分布を推定する。

好ましくは、前記分析部は、前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する位相スペクトルを用いて、テラヘルツ帯の屈折率スペクトルを推定する。

好ましくは、前記分析部は、前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する複素振幅スペクトルを用いて、テラヘルツ帯の吸収係数スペクトルを推定する。

好ましくは、前記分析装置は、前記レーザ光の偏光方向及び前記テラヘルツ波の偏光方向の組み合わせである複数の偏光配置のそれぞれについて、前記試料にレーザ光を入射させて、前記試料から出射した前記テラヘルツ波を検出し、前記分析部は、前記照射位置ごとに、複数の前記偏光配置それぞれの前記観測データを時間軸方向に並べることで合成された波形ベクトルを生成し、合成された波形ベクトルに対して統計的データ解析を行うことによって、前記ドメインの分布を推定する。

また、本発明にかかる情報処理装置は、試料の第１の面に沿ってレーザ光を走査して前記レーザ光を前記第１の面の側から前記試料に入射させ、前記試料から出射したテラヘルツ波を検出する分析装置を制御する分析装置制御部と、前記テラヘルツ波の観測された振幅の時間的な推移を示す観測データを取得する取得部と、取得された前記観測データに対して統計的データ解析を行うことによって、前記試料の、前記レーザ光の入射方向に沿った方向である深さ方向の複数の深さ位置それぞれにおける、前記レーザ光を走査した位置に対応する照射位置に関するドメインの分布を推定する分析部と、を有する。

また、本発明にかかる分析方法は、試料の第１の面に沿ってレーザ光を走査して前記レーザ光を前記第１の面の側から前記試料に入射させ、前記試料から出射したテラヘルツ波を検出し、前記テラヘルツ波の観測された振幅の時間的な推移を示す観測データを取得し、取得された前記観測データに対して統計的データ解析を行うことによって、前記試料の、前記レーザ光の入射方向に沿った方向である深さ方向の複数の深さ位置それぞれにおける、前記レーザ光を走査した位置に対応する照射位置に関するドメインの分布を推定する。

好ましくは、前記試料の前記複数の深さ位置それぞれで発生し前記試料内を伝搬した後に前記試料の外に出射された前記テラヘルツ波に関する前記観測データを取得し、取得された前記観測データに対して統計的データ解析を行うことによって、取得された前記観測データに対応する前記テラヘルツ波の波形を推定する。

好ましくは、前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する位相スペクトルを用いて、推定された前記ドメインの分布に対応する前記深さ位置を推定する。

好ましくは、取得された前記観測データに対して因子分析を行うことによって、前記ドメインの分布を推定する。

好ましくは、前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する位相スペクトルを用いて、テラヘルツ帯の屈折率スペクトルを推定する。

好ましくは、前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する複素振幅スペクトルを用いて、テラヘルツ帯の吸収係数スペクトルを推定する。

好ましくは、前記レーザ光の偏光方向及び前記テラヘルツ波の偏光方向の組み合わせである複数の偏光配置のそれぞれについて、前記試料にレーザ光を入射させて、前記試料から出射した前記テラヘルツ波を検出し、前記照射位置ごとに、複数の前記偏光配置それぞれの前記観測データを時間軸方向に並べることで合成された波形ベクトルを生成し、合成された波形ベクトルに対して統計的データ解析を行うことによって、前記ドメインの分布を推定する。

また、本発明にかかるプログラムは、試料の第１の面に沿ってレーザ光を走査して前記レーザ光を前記第１の面の側から前記試料に入射させ、前記試料から出射したテラヘルツ波を検出する分析装置を制御するステップと、前記テラヘルツ波の観測された振幅の時間的な推移を示す観測データを取得するステップと、取得された前記観測データに対して統計的データ解析を行うことによって、前記試料の、前記レーザ光の入射方向に沿った方向である深さ方向の複数の深さ位置それぞれにおける、前記レーザ光を走査した位置に対応する照射位置に関するドメインの分布を推定するステップと、をコンピュータに実行させる。

本発明によれば、試料の結晶内部の３次元的な性質を効率的に可視化することが可能な分析システム、情報処理装置、分析方法及びプログラムを提供できる。

実施の形態１にかかる分析システムの構成を示す図である。 実施の形態１にかかる試料保持装置を示す図である。 実施の形態１にかかる試料保持装置の動作を説明するための図である。 実施の形態１にかかる情報処理装置の構成を示す図である。 実施の形態１にかかる分析システムによって実行される分析方法を示すフローチャートである。 実施の形態１にかかる、試料で発生するテラヘルツ波形を説明するための図である。 実施の形態１にかかる、試料で発生するテラヘルツ波形を説明するための図である。 実施の形態１にかかる、試料で発生するテラヘルツ波形を説明するための図である。 実施の形態１において、観測されるテラヘルツ波の波形を離散的に示す図である。 試料のある照射位置にレーザパルスを照射したときに発生したテラヘルツ波形を例示する図である。 図１０に示したテラヘルツ波形のスペクトル（周波数スペクトル）を例示する図である。 実施の形態１にかかる因子分析の結果を例示する図である。 実施の形態１にかかる因子分析の結果を例示する図である。 実施の形態１にかかる因子分析の結果を例示する図である。 図１２に示したテラヘルツ波形に対応する位相スペクトルを例示する図である。 図１２～図１５に示した例における、最も位相差の大きい２組の因子成分の波形を示す図である。 図１２～図１５に示した例における、最も位相差の大きい２組の因子成分の位相スペクトルを示す図である。 図１２～図１５に示した例における、因子分析の結果を用いて推定されたテラヘルツ帯の屈折率スペクトルを例示する図である。 図１２～図１５に示した例における、因子分析の結果を用いて推定されたテラヘルツ帯の吸収係数スペクトルを例示する図である。 実施の形態１にかかる、信頼できる因子数の決定方法を説明するための図である。 テラヘルツ振幅の絶対値の３次元行列のヒストグラムを示す図である。 本実施の形態にかかるドメインイメージング方法と他の光学的ドメインイメージング方法との比較を示す図である。 実施の形態２にかかる、テラヘルツ波発生の偏光依存性測定の概念を説明するための図である。 実施の形態２にかかる、各偏光配置で観測されたテラヘルツ波形を例示する図である。 実施の形態２にかかる因子分析の方法を説明するための図である。 実施の形態２にかかる分析システムによって実行される分析方法を示すフローチャートである。 実施の形態３にかかる光学配置の例を示す図である。 実施の形態４にかかる方法を説明するための図である。

（実施の形態の概要）
従来の３次元的強誘電ドメイン可視化技術の問題点について述べる。なお、ここでは簡単のため強誘電体に限定して記述するが、本実施の形態で説明する手法自体は強磁性体及び焦電体などにも応用可能である。

３次元的強誘電ドメインの可視化においては、非特許文献１～４に示すような第二高調波発生が研究されている（その詳細は後述する）。しかしながら、対象物質に求められる位相整合条件が厳しく、また、表面の凸凹が１００ｎｍ（近赤外レーザの第二高調波の波長の１／１０程度）以下であることが必要であるため、対象物の光学研磨が必要であった。さらに、対象物質について、位相整合条件が満たされていることを担保するため屈折率スペクトル等の光学スペクトル情報（分光データ）を予め準備しておく必要がある。したがって、第二高調波発生では、試料の結晶内部の３次元的強誘電ドメインを効率的に可視化することができないおそれがある。

また、テラヘルツ波発生によって強誘電ドメイン構造を可視化する研究を、本発明者は報告している（例えば非特許文献５）。この方法は、レーザ光が結晶内部に侵入することができるならば、発生する電磁波（テラヘルツ波）の波長の約１／１０の１０μｍ程度まで表面粗さにロバストである。

本発明者は、非特許文献５において、複数の偏光方向のテラヘルツ波が発生し、かつそれぞれの偏光のテラヘルツ波に対して吸収係数が大きく異なる結晶において、検出するテラヘルツ波の偏光方向を変えることによって擬三次元的にドメイン構造を可視化できることを報告している。しかしながら、非特許文献５にかかる技術では、テラヘルツ波のピーク振幅のみを測定しており、テラヘルツ波に対する「再吸収」の大小によって、これが小さければ深さ方向で平均した像、大きければテラヘルツ波出射側の面付近の像を見ていることを利用している。したがって、ドメイン構造の深さ方向（電磁波の進行方向）の分割した可視化を実現することは困難である。

上述したような、従来のテラヘルツ波発生による強誘電ドメインイメージングでは、レーザスポット位置をスキャンしながら、ピーク振幅のみを観測していた。これに対し、本発明者は、レーザ光の各照射位置で照射したときに観測されたテラヘルツ波を統計的に解析することにより、深さ方向で断面構造を可視化できることを着案した。着眼した点は、それぞれの深さで発生したテラヘルツ波の「和」として観測されるテラヘルツ波が結晶外に出射しているということである。そして、さらに、本発明者は、ある深さで発生し放射されるテラヘルツ波の波形はその深さで一定であることを見出した。さらに、本発明者は、屈折率などの光学スペクトルの情報がない状態でもこの手法は有効であることを見出した。

（実施の形態１）
以下、本発明を適用した具体的な実施形態について、図面を参照しながら詳細に説明する。但し、本発明が以下の実施形態に限定される訳ではない。また、説明を明確にするため、以下の記載及び図面は、適宜、簡略化されている。

図１は、実施の形態１にかかる分析システム１の構成を示す図である。分析システム１は、分析装置１０及び情報処理装置１００を有する。分析装置１０は、試料９０の入射面９０ａ（第１の面）に沿ってレーザパルスを走査してレーザパルス（レーザ光）を入射面の側から試料に入射させる。そして、分析装置１０は、試料９０から出射したテラヘルツ波を検出する。詳しくは後述する。

分析装置１０は、例えば、レーザ走査テラヘルツ放射顕微鏡（Laser scanning THz emission microscopy（ＬＴＥＭ））である。分析装置１０は、レーザ光源１２と、ビームスプリッタ１４と、オプティカルチョッパ１６と、反射型可変ＮＤ（Neutral Density）フィルタ１８と、１／２波長板２０と、近赤外対物レンズ２２と、試料保持装置５０とを有する。レーザ光源１２は、フェムト秒レーザパルス２（レーザパルス）を射出する。

また、分析装置１０は、軸外し放物面鏡２４，３０と、黒色低密度ポリエチレンフィルム２６と、ワイヤーグリッド偏光子２８Ａ，２８Ｂ（wire-grid polarizer：ＷＧ）とを有する。黒色低密度ポリエチレンフィルム２６（Black LDPE（Low Density Polyethylene））の厚さは、例えば５μｍである。また、分析装置１０は、光伝導スイッチ３２と、光学遅延ステージ４０と、レンズ４４と、電流電圧変換器４６と、ロックイン検出器４８とを有する。

試料保持装置５０は、結晶である試料９０を保持する。試料保持装置５０に保持された試料９０にフェムト秒レーザパルス２が照射されると、試料９０からテラヘルツ波が出射する。ここで、本実施の形態では、右手系の座標系を考える。具体的には、試料９０にフェムト秒レーザパルス２が入射する方向（光軸方向）を、Ｚ軸方向とする。そして、Ｚ軸正方向を見て左方向をＸ軸方向とし、上方向をＹ軸方向とする。そして、試料９０の光軸方向の位置（試料９０の深さ方向の位置）をｚで示す。また、試料９０のＺ軸に垂直な面内（以下、単に面内又はｘｙ面内と称する）のＸ軸方向の座標をｘとし、試料９０の面内のＹ軸方向の座標をｙとする。また、試料９０の厚さ（光軸方向（深さ方向）の寸法）をｄとする。

分析装置１０において、レーザ光源１２から射出されたフェムト秒レーザパルス２は、ビームスプリッタ１４で、ポンプパルス４とトリガパルス６（プローブパルス）とに分割される。ここで、レーザ光源１２は、例えば、モードロックＴｉ：Ｓａｐｐｉｒｅレーザ（Mode-locked Ti:sapphire laser）である。レーザ光源１２から射出されたフェムト秒レーザパルス２のパルス幅は、例えば、１３０ｆｓである。なお、以降、フェムト秒レーザパルス２を単にレーザパルスと称することがある。

ポンプパルス４は、オプティカルチョッパ１６によって、例えば４．５ｋＨｚに強度変調される。その強度変調された周期信号は、電気回路によって参照信号としてロックイン検出器４８に入力される。また、ポンプパルス４は、反射型可変ＮＤフィルタ１８と、１／２波長板２０と、近赤外対物レンズ２２とを通過して、試料保持装置５０に保持された試料９０に照射される。これにより、入射面９０ａの側から試料９０にフェムト秒レーザパルス２（ポンプパルス４）が入射する。これにより、試料９０の、入射面９０ａとは反対側の出射面９０ｂ（第２の面）から、テラヘルツ波８０（テラヘルツ電磁波）が出射する。なお、試料保持装置５０の構成については後述する。

試料９０から出射したテラヘルツ波８０は、軸外し放物面鏡２４で反射し、黒色低密度ポリエチレンフィルム２６とワイヤーグリッド偏光子２８Ａ（ＷＧ１），２８Ｂ（ＷＧ２）とを通過し、軸外し放物面鏡３０で反射して、光伝導スイッチ３２によって受光される。なお、光伝導スイッチ３２（光伝導アンテナ）は、例えば、低温成長ＧａＡｓ光伝導スイッチである。光伝導スイッチ３２は、テラヘルツ波８０を検出するために用いられる。

トリガパルス６は、レンズ４４を通過して光伝導スイッチ３２に入射する。ここで、トリガパルス６の光路上には、光学遅延ステージ４０が設けられている。光学遅延ステージ４０は、情報処理装置１００によって制御される。光学遅延ステージ４０によってポンプパルス４とトリガパルス６との間に光学遅延を与えることにより、その光学遅延に応じた時刻のテラヘルツ波８０の電場振幅を検出することができる。可動鏡４２を移動させることで光路長が変化するので、トリガパルス６が光伝導スイッチ３２に到達する時間を変化させることができる。これにより、トリガパルス６の到達タイミングごとにサンプリングを行うことで、光伝導スイッチ３２において、テラヘルツ波８０の時間波形が得られる。光伝導スイッチ３２は電流出力であり、光伝導スイッチ３２によって得られたテラヘルツ波８０の時間波形は、電流電圧変換器４６で電圧に変換されて増幅される。ロックイン検出器４８（ロックインアンプ）は、オプティカルチョッパ１６からの参照信号に基づいて、テラヘルツ波８０の時間波形を抽出（検出）する。

情報処理装置１００は、例えばコンピュータである。情報処理装置１００は、少なくとも、試料保持装置５０及び光学遅延ステージ４０の動作を制御する。また、情報処理装置１００は、１／２波長板２０及びワイヤーグリッド偏光子２８を制御してもよい。情報処理装置１００は、ロックインアンプから、テラヘルツ波に関する情報を受信する。そして、情報処理装置１００は、テラヘルツ波に関する情報を用いて、試料９０の結晶内部の３次元的な性質を可視化するための処理を行う。情報処理装置１００の構成については後述する。なお、「結晶内部の３次元的な性質」とは、結晶の対称性を特徴づける要素であって、例えば分極（電気分極又は磁気分極）の向き及び大きさである。また、各ドメインそれぞれが有する性質は、そのドメイン内部で実質的に同じであり、各ドメインは、他のドメインとは異なる性質を有し得る。

図２は、実施の形態１にかかる試料保持装置５０を示す図である。また、図３は、実施の形態１にかかる試料保持装置５０の動作を説明するための図である。試料保持装置５０は、試料保持部５２と、可動ステージ５４，５６とを有する。試料保持部５２は、試料９０を保持する。具体的には、試料保持部５２は、試料９０の入射面９０ａをＺ軸の負方向に、出射面９０ｂをＺ軸の正方向に向けるようにして、試料９０を保持する。これにより、試料９０の入射面９０ａにフェムト秒レーザパルス２（ポンプパルス４）が入射され、出射面９０ｂからテラヘルツ波８０が出射する。

可動ステージ５４は、矢印Ａｙで示すように、試料保持部５２をＹ軸に沿った方向に移動させるように構成されている。可動ステージ５６は、矢印Ａｘで示すように、試料保持部５２をＸ軸に沿った方向に移動させるように構成されている。可動ステージ５４，５６が動作することによって、図３の矢印で示すように、フェムト秒レーザパルス２が、試料９０の入射面９０ａを、ＸＹ平面に沿ってラスタースキャンする。そして、ｋをｘｙ面内でのレーザ照射位置のインデックスとすると、分析装置１０は、ｘｙ面内の各位置ｋそれぞれに対してフェムト秒レーザパルス２を照射して、発生したテラヘルツ波８０を観測する。なお、照射位置ｋ（レーザスポット位置ｋ）と、ｘｙ面における座標（ｘ，ｙ）とは、予め対応付けられている。そして、情報処理装置１００は、位置ｋと座標（ｘ，ｙ）との対応を示すデータを、予め記憶している。

ここで、位置ｋの総数をＫとする。したがって、ｋは、１以上Ｋ以下の整数である。フェムト秒レーザパルス２が試料９０の入射面９０ａをラスタースキャンすることによって、Ｋ個のテラヘルツ波形を測定できる。したがって、Ｋは測定される波形数を示す。波形数Ｋは、得たい面内分解能のピクセル数に対応する。例えば、直径１０μｍのレーザスポット径で、１ｍｍ×１ｍｍのエリアをイメージングする場合、１ｍｍ／１０μｍ＝１００であるので、Ｋ＝１００×１００＝１００００個のテラヘルツ波形を観測できる。

図４は、実施の形態１にかかる情報処理装置１００の構成を示す図である。また、図５は、実施の形態１にかかる分析システム１によって実行される分析方法を示すフローチャートである。なお、図５については後述する。情報処理装置１００は、主要なハードウェア構成として、制御部１０２と、記憶部１０４と、通信部１０６と、インタフェース部１０８（ＩＦ；Interface）とを有する。制御部１０２、記憶部１０４、通信部１０６及びインタフェース部１０８は、データバスなどを介して相互に接続されている。

制御部１０２は、例えばＣＰＵ（Central Processing Unit）等のプロセッサである。制御部１０２は、制御処理及び演算処理等を行う演算装置としての機能を有する。記憶部１０４は、例えばメモリ又はハードディスク等の記憶デバイスである。記憶部１０４は、例えばＲＯＭ（Read Only Memory）又はＲＡＭ（Random Access Memory）等である。記憶部１０４は、制御部１０２によって実行される制御プログラム及び演算プログラム等を記憶するための機能を有する。また、記憶部１０４は、処理データ等を一時的に記憶するための機能を有する。記憶部１０４は、データベースを含み得る。

通信部１０６は、他の装置と有線又は無線のネットワーク等を介して通信を行うために必要な処理を行う。通信部１０６は、通信ポート、ルータ、ファイアウォール等を含み得る。インタフェース部１０８は、例えばユーザインタフェース（ＵＩ；User Interface）である。インタフェース部１０８は、キーボード、タッチパネル又はマウス等の入力装置と、ディスプレイ又はスピーカ等の出力装置とを有する。インタフェース部１０８は、ユーザ（オペレータ等）によるデータの入力の操作を受け付け、ユーザに対して情報を出力する。

また、情報処理装置１００は、その機能を実現する構成要素として、分析装置制御部１１０、取得部１２０、分析部１３０及び分析結果出力部１４０を有する。なお、各構成要素は、例えば、制御部１０２の制御によって、プログラムを実行させることによって実現できる。より具体的には、各構成要素は、記憶部１０４に格納されたプログラムを、制御部１０２が実行することによって実現され得る。また、必要なプログラムを任意の不揮発性記録媒体に記録しておき、必要に応じてインストールすることで、各構成要素を実現するようにしてもよい。また、各構成要素は、プログラムによるソフトウェアで実現することに限ることなく、ハードウェア、ファームウェア、及びソフトウェアのうちのいずれかの組み合わせ等により実現してもよい。また、各構成要素は、例えばＦＰＧＡ（field-programmable gate array）又はマイコン等の、ユーザがプログラミング可能な集積回路を用いて実現してもよい。この場合、この集積回路を用いて、上記の各構成要素から構成されるプログラムを実現してもよい。

分析装置制御部１１０は、分析装置１０を制御して、試料９０の位置ｋにレーザパルスを照射させる。これにより、分析装置１０は、試料９０から発生するテラヘルツ波８０を観測する。具体的には、分析装置制御部１１０は、レーザ光源１２及び試料保持装置５０を制御して、位置ｋに、フェムト秒レーザパルス２を照射させる。また、分析装置制御部１１０は、光学遅延ステージ４０を制御して、トリガパルス６が光伝導スイッチ３２に到達する時間を制御する。これにより、位置ｋで発生したテラヘルツ波８０の時間波形が得られる。そして、分析装置制御部１１０は、試料保持装置５０を制御して、次の位置ｋ＋１にフェムト秒レーザパルス２が照射されるように、試料９０を移動させる。このような動作を繰り返すことによって、試料９０にフェムト秒レーザパルス２が走査される。これにより、全ての位置ｋ（ｋ＝１～Ｋ）について、試料９０から出射したテラヘルツ波８０が検出される。

取得部１２０は、上述した方法によって観測されたテラヘルツ波８０を示す観測データを取得する。分析部１３０は、取得された観測データを用いて後述するような分析を行う。分析結果出力部１４０は、分析結果を出力する。取得部１２０、分析部１３０及び分析結果出力部１４０の処理については、図５を用いて後述する。

以下、観測されたテラヘルツ波形を用いて試料９０の３次元的な性質（分極ドメイン）を推定する方法について説明する。
［テラヘルツ波について］
図６～図８は、実施の形態１にかかる、試料９０で発生するテラヘルツ波形を説明するための図である。図６～図８では、試料９０の深さ方向（Ｚ軸方向；光軸方向）に着目している。図６に示すように、ポンプパルス４を入射面９０ａの側から試料９０に照射したときに、ある深さ方向の点（ｚ，ｚ＋ｄｚ）で、テラヘルツ波８２が発生したとする。なお、ｄｚは、深さ方向の微小区間である。

このとき、テラヘルツ波８２の波形は、テラヘルツ波８２の試料９０への吸収、及び、入射波（ポンプパルス４）と出射波（テラヘルツ波８２）との間の位相整合条件によって、変形及び時間遅延を生じる。また、検出器（光伝導スイッチ３２，ロックイン検出器４８）で観測されるテラヘルツ波８２の複素スペクトルＥ（ω，ｚ）は、以下の式１で表される。なお、この複素スペクトルＥ（ω，ｚ）は、深さ依存性のない要素に影響される装置関数Ｈ（ω）を乗算することによって得られる。深さ依存性のない要素とは、例えば、テラヘルツ波８２の出射面９０ｂでの透過率、テラヘルツ光路での回折、検出器の応答関数、及びレーザパルス幅がある。

・・・（１）

ここで、ωは、テラヘルツ波８２の角周波数であり、テラヘルツ波８２の周波数をｆとすると、ω＝２πｆである。Ｅ（ω，ｚ）は、深さｚで発生して結晶（試料９０）を伝搬し結晶外に放射される、角周波数ωのテラヘルツ波８２の周波数空間での複素振幅である。χ^（２）（ｚ，ω）は、深さｚで発生したテラヘルツ波８２の発生機構における二次の非線形光学係数である。ｎ_ｇは、レーザパルスに対する群屈折率を示す。α_{ｌａｓｅｒ}は、レーザパルスに対する吸収係数である。ｎ（ω）は、テラヘルツ帯の屈折率スペクトルを示す。α（ω）は、テラヘルツ帯の吸収係数スペクトルを示す。

また、Ｈ（ω）は、テラヘルツ波８２の振幅測定の際の測定系の応答関数である。ここで、装置関数Ｈ（ω）は、複素スペクトルとなる。なお、上述したように、Ｈ（ω）には、例えば、テラヘルツ波８２の回折及び透過率、光伝導スイッチ３２の光キャリア寿命、レーザパルス幅、及び、レーザ光強度などが影響する。Ｈ（ω）を逆フーリエ変換することにより、時間領域の装置関数Ｈ（ｔ）が得られる。ここで、装置関数Ｈ（ω）を実験的に決定することは困難である。

そして、式１を逆フーリエ変換（ＩＦＦＴ：Inverse Fast Fourier Transform）することで得られる以下の式２が、試料９０の深さ方向の各点（ｚ，ｚ＋ｄｚ）で発生したテラヘルツ波８２を実際に分析装置１０（テラヘルツ波形検出系）によって観測して得られる観測結果に対応する。

・・・（２）

ここで、Ｅ（ｔ，ｚ）は、深さｚで発生して結晶を伝搬し結晶外に放射されたテラヘルツ波８２の時間領域での電場振幅波形を示す。なお、実際に光伝導スイッチ３２で観測される波形は、式２で表される各深さｚで発生したテラヘルツ波８２の波形についてｚ＝０～ｄで和をとり、さらに時間領域の装置関数Ｈ（ｔ）を畳み込み積分したものに対応する。

また、式２で表される波形の振幅は、二次の非線形光学係数χ^（２）（ｚ）や分極値Ｐ（ｚ）に比例する。なお、Ｐ（ｘ，ｙ，ｚ）は、結晶内の位置（ｘ，ｙ，ｚ）における強誘電分極を示す。なお、フーリエ変換の基本原理として、波形の時間遅れは、周波数空間では、位相の線形関数の積として表される。つまり、時間波形でのΔｔ（＞０）の時間差（時間遅れ）は、ｅｘｐ（－ｉΔｔω）の複素スペクトルへの乗算に対応する。なお、Δｔは、テラヘルツ波形測定での時間差を示す。

ここで、上述したように、実際に観測されるテラヘルツ波８０の波形は、試料９０のそれぞれの深さで発生したテラヘルツ波８２の和（線形和）として表される。つまり、試料９０のそれぞれの深さで発生したテラヘルツ波８２が重ね合わされたテラヘルツ波８０の波形が結晶外に出射している。以下、図７及び図８を用いて説明する。

図７は、均質な非線形光学結晶である試料９０において、フェムト秒レーザパルス２（ポンプパルス４）によって発生したテラヘルツ波８２の伝播を概念的に示す。時刻ｔ１において、フェムト秒レーザパルス２が深さｚ１に到達すると、深さｚ１において、テラヘルツ波８２－１が発生する。そして、テラヘルツ波８２－１は、Ｚ軸の正方向に伝搬する。時刻ｔ２において、フェムト秒レーザパルス２が深さｚ２に到達すると、深さｚ２において、テラヘルツ波８２－２が発生する。そして、テラヘルツ波８２－２は、Ｚ軸の正方向に伝搬する。時刻ｔ３において、フェムト秒レーザパルス２が深さｚ３に到達すると、深さｚ３において、テラヘルツ波８２－３が発生する。そして、テラヘルツ波８２－３は、Ｚ軸の正方向に伝搬する。時刻ｔ４において、フェムト秒レーザパルス２が深さｚ４に到達すると、深さｚ４において、テラヘルツ波８２－４が発生する。そして、テラヘルツ波８２－４は、Ｚ軸の正方向に伝搬する。時刻ｔ５において、フェムト秒レーザパルス２が深さｚ５に到達すると、深さｚ５において、テラヘルツ波８２－５が発生する。そして、テラヘルツ波８２－５は、Ｚ軸の正方向に伝搬する。

ここで、深さｚ１～ｚ５で同位相でそれぞれ発生したテラヘルツ波８２－１～８２－５は、テラヘルツ帯の屈折率に対応する速さで伝搬する。なお、試料９０の内部では、テラヘルツ波の速度は、フェムト秒レーザパルス２の速度よりも遅い。したがって、図７に示すように、先に発生したテラヘルツ波８２が、後で発生したテラヘルツ波８２の後ろを進むこととなる。つまり、試料９０の入射面９０ａに近い深さで発生したテラヘルツ波８２は、試料９０の出射面９０ｂに近い深さで発生したテラヘルツ波８２よりも、遅延することとなる。そして、深さｚ１～ｚ５で発生したテラヘルツ波８２－１～８２－５の総和が、テラヘルツ波８０として観測される。
ここで、試料９０の内部で、テラヘルツ波の位相伝播速度がフェムト秒レーザパルス２の強度包括線の伝播速度よりも遅い場合について述べた。しかしながら、テラヘルツ波の位相伝播速度がフェムト秒レーザパルス２の強度包括線の伝播速度よりも速い場合についても分析への影響はない。この場合については後述する。

図８は、光軸方向に分極ドメイン構造を有する強誘電体である試料９０において、フェムト秒レーザパルス２（ポンプパルス４）によって発生したテラヘルツ波８２の伝播を概念的に示す。時刻ｔ１において、フェムト秒レーザパルス２が深さｚ１に到達すると、深さｚ１において、テラヘルツ波８２－１が発生する。そして、テラヘルツ波８２－１は、Ｚ軸の正方向に伝搬する。時刻ｔ２において、フェムト秒レーザパルス２が深さｚ２に到達すると、深さｚ２において、テラヘルツ波８２－２が発生する。そして、テラヘルツ波８２－２は、Ｚ軸の正方向に伝搬する。時刻ｔ３において、フェムト秒レーザパルス２が深さｚ３に到達すると、深さｚ３において、テラヘルツ波８２－３が発生する。そして、テラヘルツ波８２－３は、Ｚ軸の正方向に伝搬する。時刻ｔ４において、フェムト秒レーザパルス２が深さｚ４に到達すると、深さｚ４において、テラヘルツ波８２－４が発生する。そして、テラヘルツ波８２－４は、Ｚ軸の正方向に伝搬する。時刻ｔ５において、フェムト秒レーザパルス２が深さｚ５に到達すると、深さｚ５において、テラヘルツ波８２－５が発生する。そして、テラヘルツ波８２－５は、Ｚ軸の正方向に伝搬する。

図７の例と同様に、深さｚ１～ｚ５で同位相でそれぞれ発生したテラヘルツ波８２－１～８２－５は、テラヘルツ帯の屈折率に対応する速さで伝搬する。ここで、各深さｚで発生するテラヘルツ波８２の振幅は、その深さｚにおける分極Ｐに依存（比例）する。そして、図８の例では、深さｚ１～ｚ５における分極Ｐ_１～Ｐ_５は、それぞれ、＋Ｐ，－Ｐ，＋Ｐ，－Ｐ，＋Ｐである。したがって、深さｚ２及びｚ４で発生したテラヘルツ波８２－２，８２－４では、符号が反転している。したがって、図８の例における、深さｚ１～ｚ５で発生したテラヘルツ波８２－１～８２－５の総和として観測されるテラヘルツ波８０の波形は、図７の場合と異なる。つまり、観測されるテラヘルツ波８０の波形は、試料９０の各深さ位置の分極構造に依存することとなる。

なお、試料９０の入射面９０ａで発生したテラヘルツ波８２－１のピーク位置と、試料９０の出射面９０ｂで発生したテラヘルツ波８２－５のピーク位置との、波形測定における時間差をΔｔとすると、Δｔは、概ね、ｃΔｎ／ｄに一致する。ここで、ｃは光速を示す。ｃは、真空中の光速であってもよい。また、Δｎは、テラヘルツ波の屈折率とレーザパルスの群屈折率との差を示す。なお、テラヘルツ波の屈折率を正に、レーザパルスの群屈折率を負にとる。

図９は、実施の形態１において、観測されるテラヘルツ波８０の波形を離散的に示す図である。実際に観測されるテラヘルツ波８０は、後述する式３で示すように、離散的なＭ点の成分からなるベクトルとして表される。つまり、テラヘルツ波８０は、Ｍ個の点でサンプリングされたデータとして観測される。ここで、各測定時刻の間隔をΔｔ_ｅｘｐ、テラヘルツ波８０の始点の時刻及び最終点の時刻をそれぞれｔ_ｍｉｎ，ｔ_ｍａｘとすると、ｔ_ｍａｘ－ｔ_ｍｉｎ＝ＭΔｔ_ｅｘｐである。また、離散化されたテラヘルツ波形での時刻インデックスをｍ（ｍは１以上Ｍ以下の整数）とする。このとき、ｔ_ｍｉｎはｔ_１であり、ｔ_ｍａｘはｔ_Ｍである。なお、Ｍは、テラヘルツ波形測定での総点数である。また、Δｔ_ｅｘｐは、Δｔの実験での設定値である。Δｔ_ｅｘｐは、例えば、０．０５～０．１ｐｓに設定される。なお、フーリエ変換の定理から、Δｔ_ｅｘｐを、測定系の最大周波数（例えば３．０ＴＨｚ）の２倍の逆数（１／６．０ｐｓ＝０．１６６ｐｓ）以下に設定すれば、テラヘルツ波８０の信号を良好に再現できる。

ここで、ｘｙ面内位置ｋにレーザパルスを照射することで発生したテラヘルツ波８０をベクトル（波形ベクトル）として表現すると、以下の式３で表される。ここで、Ｅ_ｋｍは、面内位置ｋで発生したテラヘルツ波８０の、時刻ｔ_ｍでの振幅値を示す。式３で示すように、ベクトルｅ_ｋは、サンプリング数Ｍ個の成分で表される。

・・・（３）

また、ｘｙ面内の全ての位置（ｋ＝１～Ｋ）について、テラヘルツ波８０の振幅を行列Ｅ_ＴＨｚとして表現すると、以下の式４で表される。なお、式４の行列（「振幅行列」と称する）は、Ｋ×Ｍの２次元行列である。また、式４の行列は、全ての照射点（ｋ＝１～Ｋ）について観測されるテラヘルツ波形を表した行列となる。なお、式中において、太字は行列（又はベクトル）を示す（他の式についても同様）。

・・・（４）

ここで、試料９０の深さ方向（光軸方向；Ｚ軸方向）の位置のインデックスをｌ（ｌは１以上Ｌ以下の整数）とする。そして、深さ位置ｌ（深さｚ_ｌ（０≦ｚ_ｌ≦ｄ））で発生したテラヘルツ波８２の波形ベクトルをｖ_ｌとする。なお、波形ベクトルｖ_ｌは、基底ベクトルであるとする。波形ベクトルｖ_ｌは、式３と同様に、Ｍ個の成分で表され得る。なお、後述するように、観測したテラヘルツ波８０から波形ベクトルｖ_ｌの各成分値が直接測定できるわけではない。また、波形ベクトルｖ_ｌは、ｘｙ面内位置ｋとは独立している。つまり、波形ベクトルｖ_ｌは、ｘｙ面内位置ｋに依存せず、深さ位置ｌのみに依存する。また、レーザパルスの各照射位置ｋのときに観測されるテラヘルツ波８０の波形は、この波形ベクトルｖ_ｌで示される波形をｌ＝１～Ｌについて重ね合わせたものに対応する。そして、図８を用いて説明したように、レーザパルスの照射位置ｋのときに深さ位置ｌで発生したテラヘルツ波８２の振幅は、試料９０内の点（ｋ，ｌ）での分極値Ｐ_ｋｌに比例する。

以下の式５は、各点（ｋ，ｌ）の分極値Ｐ_ｋｌを行列で表したものである。なお、行列Ｐ_ｋｌは、Ｋ×Ｌの２次元行列である。なお、以下、行列Ｐ_ｋｌを、分極行列と称することがある。

・・・（５）

例えば、ベクトルＰ_ｋ１（式５の１列目の列ベクトルに対応）は、試料９０の深さ方向の１層目（ｌ＝１）のｘｙ面内におけるドメイン構造を表す分極のベクトルである。試料９０が一様な非線形光学感受率を持つ非線形光学結晶である場合は、Ｐ_ｋ１が一定のもの（つまりＰ_ｋ１の各成分が互いに同じ）としてモデル化することができる。

ここで、深さ方向の位置の分割数をＬとする。上述したように、レーザパルスの照射位置ｋのときに深さ位置ｌで発生したテラヘルツ波８２の振幅は、その点（ｋ，ｌ）での分極値Ｐ_ｋｌに比例する。そして、各レーザ照射位置ｋについて観測されるテラヘルツ波８０は、分極値Ｐ_ｋｌを振幅とし、波形ベクトルｖ_ｌで示される波形をｌ＝１～Ｌについて重ね合わせたものに対応する。したがって、各レーザ照射位置ｋについて観測されるテラヘルツ波８０は、以下の式６で表される。

・・・（６）

なお、式６を行列のまま表現すると、以下の式７のように表される。なお、行列Ｅ_{ＴＨｚ，ｋ}は、式４と同様に、Ｋ×Ｍの２次元行列である。

・・・（７）

式６（及び式７）で示すように、あるスポット位置ｋにレーザパルスを照射したときに観測されるテラヘルツ波形（行列Ｅ_{ＴＨｚ，ｋ}）は、各深さ位置ｌからの寄与の線形和で表される。そして、実施の形態１にかかる３次元的イメージングでは、分極Ｐ_ｋｌをスポット位置ｋ及び深さ位置ｌの双方について推定する。

しかしながら、式６，式７に示したｖ_ｌを実験的に求めることは、以下に説明するように、極めて困難である。ｖ_ｌを完全に解析的に決定するためには、テラヘルツ帯の屈折率スペクトルｎ（ω）及び吸収係数スペクトルα（ω）、レーザ光に対する群屈折率ｎ_ｇ及び吸収係数α_{ｌａｓｅｒ}、さらに、テラヘルツ波形計測の装置関数Ｈ（ω）が必要である。とくにＨ（ω）の決定は非常に困難であり、Ｈ（ω）を推定できたとしてもその測定誤差がＨ（ω）に影響を与える。また、大きさ１ｍｍ以上の結晶合成が難しい対象物については、テラヘルツ帯のスペクトル（ｎ（ω），α（ω））を測定することが難しい。

そこで、本実施の形態にかかる分析部１３０は、全ての照射位置ｋ（ｋ＝１～Ｋ）について実験的に得られた振幅行列Ｅ_ＴＨｚ（式４）に対して統計データ解析を行うことにより、分極行列Ｐ_ｋｌを推定する。具体的には、分析部１３０は、統計データ解析手法の「因子分析（Factor analysis）」を用いることで、分極行列Ｐ_ｋｌと波形ベクトルｖ_ｌ（基底ベクトル）を一意に決定することができる。さらに、その推定誤差も得られる。ここで、分極行列Ｐ_ｋｌの各成分Ｐ_ｋｌは、試料９０内の点（ｋ，ｌ）での振幅に対応する。また、分極行列Ｐ_ｋｌは、推定された分極ドメインの分布に対応する。

なお、本実施の形態に関して行われた実証実験では、試料９０として、有機強誘電体クロコン酸結晶を用いた。また、試料９０の厚さをｄ＝３２０μｍとした。また、試料９０のＸＹ平面に沿った方向の寸法については、Ｘ軸方向の寸法を４００μｍとし、Ｙ軸方向の寸法を６００μｍとした。

図１０は、試料９０のある照射位置ｋにレーザパルスを照射したときに発生したテラヘルツ波形を例示する図である。また、図１１は、図１０に示したテラヘルツ波形のスペクトル（周波数スペクトル）を例示する図である。なお、式３及び図９を用いて上述したように、実際には、図１０に例示したテラヘルツ波形に対応するＭ個の離散的なデータが観測される。この観測されたデータ（観測データ）は、振幅の時間的な推移を示す。

［統計データ解析手法（因子分析）］
因子分析は、例えばＰｙｔｈｏｎのようなデータ解析プログラムで実行され得る。したがって、情報処理装置１００の制御部１０２は、Ｐｙｔｈｏｎ等のデータ解析プログラムを実行し、分析部１３０は、Ｐｙｔｈｏｎ等のデータ解析プログラムによって実現され得る。また、実施の形態１における因子分析では、波形ベクトルｖ_ｌを共通因子とし、分極行列Ｐ_ｋｌを因子負荷行列とする。この場合、分極行列Ｐ_ｋｌの各成分が、因子負荷量に対応する。また、観測データに対応する振幅行列Ｅ_{ＴＨｚ，ｋ}を観測変数とする。そして、波形ベクトルｖ_ｌ及び振幅行列Ｅ_{ＴＨｚ，ｋ}が離散的なＭ個の成分からなるベクトル（式４）で表されることから、全ての位置ｋについて、式６に対応するＫ×Ｍ個の方程式を最適に満たす分極行列Ｐ_ｋｌ（因子負荷行列）及び波形ベクトルｖ_ｌ（共通因子）を求める。より具体的には、全ての位置ｋについて、式６に対応するＫ×Ｍ個の方程式を最適に満たす分極行列Ｐ_ｋｌ（因子負荷行列）の各成分及び波形ベクトルｖ_ｌ（共通因子）の各成分を求める。また、因子分析において、深さ位置の因子数をＦとする。Ｆは、例えば、試行錯誤によって決定されてもよい。なお、Ｆ≦Ｌである。また、因子成分のインデックスをｂ（ｂは１以上Ｆ以下の整数）とする。

図１２～図１４は、実施の形態１にかかる因子分析の結果を例示する図である。図１２～図１４は、深さ位置の因子数Ｆを７つ（Ｆ＝７）とした場合の因子分析の結果を例示している。なお、図１２～図１３は、７つの因子成分（ｂ＝１～７）のそれぞれに関する波形ベクトルｖ_ｌ（ｌ＝１～７）に対応する。但し、ｌの値とｂの値とは一致するわけではないことに留意されたい。例えば、深さ位置ｌ＝１で発生した波形ベクトルｖ_１は、因子成分ｂ＝１（第１因子成分）の波形ベクトルに対応するわけではない。

図１２は、７つの因子成分のそれぞれについて、因子分析によって抽出された波形ベクトルｖ_ｌ（共通因子）のＭ個の点（成分）でサンプリングされたデータをプロットして得られたテラヘルツ波形を示す。また、図１３は、図１２に示すテラヘルツ波形それぞれのスペクトルを示す。図１２及び図１３において、７つの因子成分のそれぞれについての波形が、互いに異なる線種で示されている。例えば、第１成分（ｂ＝１）は、細い実線で示されている。なお、図１２及び図１３に示す細い二点鎖線は、ｘｙ面内のイメージ全体での平均のテラヘルツ波形を示す。なお、図１２及び図１３では、それぞれの波形及びスペクトルが見やすくなるように、各波形及び各スペクトルが縦方向にオフセット調整されて示されている。したがって、実際には、各因子成分の波形及びスペクトルの左右両端における振幅は、略０である。

図１２及び図１３から、７つの因子成分のそれぞれについて、異なる波形構造及びスペクトルであることが分かる。また、７つの因子成分のそれぞれの間には時間遅延が表れている。例えば、図１２において、第１成分（ｂ＝１）と第７成分（ｂ＝７）とでは、ピーク位置が異なっている。

また、図１４は、７つの因子成分のそれぞれに対応する深さにおける、ｘｙ面内における分極ドメイン（振幅）の分布を示す。図１４は、因子分析によって抽出された分極行列Ｐ_ｋｌ（因子負荷行列）の各成分をプロットすることで得られた分極Ｐ（ｘ，ｙ，ｚ）のドメインの分布を示している。これにより、図１４で示すように、試料９０の３次元的な分極（振幅）のイメージングを実現することができる。

図１４に示した分極ドメインの分布において、色の濃さによって、分極値（振幅）の大小を示している。色が薄い位置ほど分極値は大きく、色が濃い位置ほど分極値は小さい。また、色の濃さが黒と白との中間色（グレー）の位置では、分極値は０となっている。したがって、黒色の位置では、分極値は負である。また、図１４に示した分極ドメインの分布においては、それぞれの層（因子成分）で、カラースケールを規格化している。

なお、図１４では、どの因子成分ｂがどの深さ位置ｌに対応するかが示されている。例えば、因子成分ｂ＝１（第１成分）は、深さ位置ｌ＝７に対応する。なお、その対応付けは、因子分析を行うだけでは分からない。この、どの因子成分ｂがどの深さ位置ｌに対応するかを対応付ける方法について、以下に説明する。

図１５は、図１２に示したテラヘルツ波形に対応する位相スペクトルを例示する図である。本実施の形態にかかる分析部１３０は、この位相スペクトル（位相特性）を用いて、どの因子成分ｂがどの深さ位置ｌに対応するかを対応付ける。なお、位相スペクトルは、図１２に示した因子成分ごとのテラヘルツ波形それぞれについてフーリエ変換を行うことによって得られる。図１５に示すように、位相スペクトルは、横軸を周波数とし、縦軸を位相として表される。

図７及び図８を用いて上述したように、試料９０の入射面９０ａに近い深さで発生したテラヘルツ波８２は、試料９０の出射面９０ｂに近い深さで発生したテラヘルツ波８２よりも、遅延する。したがって、分析部１３０は、各因子成分ｂの時間遅延を分析することによって、因子成分ｂと深さ位置ｌとの対応を判定する。ここで、上述したように、フーリエ変換の基本原理として、波形の時間遅れΔｔは、周波数空間では、関数ｅｘｐ（－ｉΔｔω）の積として表される。したがって、各因子成分について位相の周波数依存性を解析することによって、各因子成分についての時間遅延を抽出することができる。

具体的には、分析部１３０は、複数の因子成分ｂ（ｂ＝１～７）についての位相スペクトルそれぞれについて、線形関数でカーブフィッティング解析（回帰分析）を行う。この線形関数（位相曲線）の傾きが、時間遅延に対応する。図１５で示すように、各因子成分ｂについて、位相の周波数依存性が異なっているので、位相曲線の傾きも異なっている。分析部１３０は、位相曲線の傾きから、因子成分ｂと深さ位置ｌとの対応を判定する。

ここで、実施の形態１にかかるＰｙｔｈｏｎのフーリエ変換の計算では、時間遅延が位相曲線（θ（ｆ）曲線）の負の傾きと対応するようになっている。したがって、位相曲線の傾きがプラス側に大きい（つまり周波数が増加すると位相がより増加する傾向にある）場合に、その因子成分ｂに対応するテラヘルツ波形の時間遅延は小さい。したがって、位相曲線の傾きがプラス側に大きい場合に、対応する因子成分ｂは、出射面９０ｂの側の深さ位置ｌに対応する。逆に言うと、位相曲線の傾きがマイナス側に大きい（つまり周波数が増加すると位相がより減少する傾向にある）場合に、その因子成分ｂに対応するテラヘルツ波形の時間遅延は大きい。したがって、位相曲線の傾きがマイナス側に大きい場合に、対応する因子成分ｂは、入射面９０ａの側の深さ位置ｌに対応する。

図１５に示すように、各因子成分ｂについて、第１成分、第２成分、第４成分、第５成分、第３成分、第６成分、第７成分の順で、位相曲線の傾きが、プラス側に大きい。つまり、第１成分（ｂ＝１）の位相曲線の傾きが最もプラス側に大きく、第７成分（ｂ＝７）の位相曲線の傾きが最もプラス側に小さい。したがって、図１５の例では、分析部１３０は、因子成分のインデックスｂ＝７は深さ位置ｌ＝１に対応すると判定する。したがって、因子成分ｂ＝７に対応する分極ドメインの分布は、分極のベクトルＰ_ｋ１（式５）に対応する。

同様に、分析部１３０は、因子成分のインデックスｂ＝６は深さ位置ｌ＝２に対応すると判定する。したがって、因子成分ｂ＝６に対応する分極ドメインの分布は、分極のベクトルＰ_ｋ２に対応する。また、分析部１３０は、因子成分のインデックスｂ＝３は深さ位置ｌ＝３に対応すると判定する。したがって、因子成分ｂ＝３に対応する分極ドメインの分布は、分極のベクトルＰ_ｋ３に対応する。また、分析部１３０は、因子成分のインデックスｂ＝５は深さ位置ｌ＝４に対応すると判定する。したがって、因子成分ｂ＝５に対応する分極ドメインの分布は、分極のベクトルＰ_ｋ４に対応する。また、分析部１３０は、因子成分のインデックスｂ＝４は深さ位置ｌ＝５に対応すると判定する。したがって、因子成分ｂ＝４に対応する分極の分布は、分極のベクトルＰ_ｋ５に対応する。

また、分析部１３０は、因子成分のインデックスｂ＝２は深さ位置ｌ＝６に対応すると判定する。したがって、因子成分ｂ＝２に対応する分極ドメインの分布は、分極のベクトルＰ_ｋ６に対応する。また、分析部１３０は、因子成分のインデックスｂ＝１は深さ位置ｌ＝７に対応すると判定する。したがって、因子成分ｂ＝１に対応する分極ドメインの分布は、分極のベクトルＰ_ｋ７に対応する。
このようにして、式５に示した分極行列の各成分Ｐ_ｋｌが求められる。したがって、分析部１３０は、試料９０の３次元的な分極ドメインを可視化することができる。

なお、図１２～図１５に示した例では、因子数Ｆ＝７であるので、深さ方向の分解能は、Δｚ＝３２０／７～４５μｍ程度である（但し「～」は近似記号）。ただし、急激に分極ドメイン構造が変化している場合、より高い深さ分解能で層が分離される。この場合、因子分析の際に有効な因子数が増え、これにより深さ分解能が向上する。また、時間分解能（最大周波数）を向上することにより、深さ方向の分解能を向上できる。また、図１５に示した位相スペクトルは、イメージ全体での平均のテラヘルツ波形を、位相計算の際の参照信号としている。これにより、深さ方向のインデックスの決定の精度を向上させることができる。

なお、本実施の形態にかかる方法における深さ方向の分解能Δｚは、下記のように決めてもよい。信号雑音比が十分大きい周波数の最大（最大周波数）をｆ_ｍａｘとする。その波長λはλ＝ｃ／ｆ_ｍａｘと表される。テラヘルツ波形測定での近接点間の時間をΔｔとすると、テラヘルツ波形測定の最大位相感度Δθは、Δθ＝ｆ_ｍａｘΔｔとなる。そして、Δｚは、Δｚ＝λΔθ＝ｃΔｔで求められる。なお、実証実験では、ｆ_ｍａｘ＝３．０ＴＨｚ（λ＝１００μｍ），Δｔ＝０．０８ｐｓであり、Δｚ＝２４μｍとなる。なお、波形測定の時間幅を長くする、またはΔｔを小さくすることにより、Δｚを改善できる。

［テラヘルツ帯の屈折率スペクトルの推定］
分析部１３０は、さらに、以下に説明するように、テラヘルツ帯の屈折率スペクトルの推定を行う。上述した各因子成分ｂに対応するテラヘルツ波形は、対応する深さ位置ｌから発生したものとみなせる。したがって、各成分間の時間遅延及び結晶厚さ（ｄ）の情報から、以下の式８で示されるΔｎ（ω）を決定することができる。

・・・（８）

ここで、Δｎ（ω）は、テラヘルツ波の屈折率スペクトルｎ（ω）とレーザ群屈折率ｎ_ｇとの差を示すスペクトルを表す。また、レーザ群屈折率ｎ_ｇは、例えばレーザ波長帯の屈折率スペクトル測定方法によって、予め取得されているものとする。したがって、式８から、Δｎ（ω）を求めてレーザ群屈折率ｎ_ｇを加算することで、テラヘルツ波の屈折率スペクトルｎ（ω）を求めることができる。

ここで、Δｎ（ω）の精度は、結晶厚さｄ及び最大観測周波数に比例する。また、最も位相差の大きい２組の因子成分の波形について、それぞれに対応する深さの差は、概ねｄ×（Ｆ－１）／Ｆとなる。したがって、Δｎ（ω）とΔθ（ω）の間には、以下の式９で示される関係式が成り立つ。なお、Δθ（ω）は、角周波数ωにおける、最も位相差の大きい２組の因子成分の位相差を示す。

・・・（９）

図１６は、図１２～図１５に示した例における、最も位相差の大きい２組の因子成分の波形を示す図である。上述した因子分析の結果、最も位相差の大きい２組の因子成分は、第１因子成分及び第７因子成分である。そして、分析部１３０は、最も位相差の大きい２組の因子成分の波形の位相スペクトルを算出し、算出された位相スペクトルから、Δθ（ω）を算出する。

図１７は、図１２～図１５に示した例における、最も位相差の大きい２組の因子成分の位相スペクトルを示す図である。分析部１３０は、各周波数について、第１因子成分の位相と第７因子成分の位相との位相差Δθを算出する。例えば、周波数ｆ＝１．０ＴＨｚのとき、Δθ＝２．５πである。また、周波数ｆ＝１．５ＴＨｚのとき、Δθ＝４．０πである。また、周波数ｆ＝２．０ＴＨｚのとき、Δθ＝５．０πである。そして、ω＝２πｆであるので、分析部１３０は、式８及び式９を用いて、各角周波数ωについて、テラヘルツ波の屈折率を算出できる。そして、分析部１３０は、各周波数とそのときの屈折率との関係から、テラヘルツ波の屈折率スペクトルｎ（ω）を算出する。

図１８は、図１２～図１５に示した例における、因子分析の結果を用いて推定されたテラヘルツ帯の屈折率スペクトルを例示する図である。図１８において、因子分析の結果を用いて推定されたテラヘルツ波の屈折率スペクトルが、黒丸で示されている。また、比較のため、本実施の形態にかかる方法とは別の方法であるテラヘルツ時間領域分光によって得られたテラヘルツ波の屈折率スペクトルが、実線で示されている。なお、図１８に示したテラヘルツ波の屈折率スペクトルは、各周波数とそのときの屈折率とをプロットすることによって得られる。因子分析の結果を用いて推定されたテラヘルツ波の屈折率スペクトルは、テラヘルツ時間領域分光によって得られたテラヘルツ波の屈折率スペクトルと、１０％程度の精度で一致している。したがって、因子分析の結果を用いて、良好に、テラヘルツ波の屈折率スペクトルを推定することができる。

［テラヘルツ帯の吸収係数スペクトルの推定］
分析部１３０は、さらに、以下に説明するように、テラヘルツ帯の吸収係数スペクトルの推定を行う。分極Ｐ_ｋｌの絶対値をｋについて平均をとったものを、以下の式１０のように表す。なお、便宜上、以後、式１０を「Ｐ^－ _ｌ」と表記することがある。

・・・（１０）

Ｐ^－ _ｌは、第ｌ層（深さ位置ｌ）から発生し検出されたテラヘルツ波８２の振幅となる。位相スペクトルθ（ω）（図１５，図１７）の傾きから計算された第ｌ層の深さ、及びα_{ｌａｓｅｒ}より、テラヘルツ帯の吸収係数スペクトルの推定を行うことができる。なお、α_{ｌａｓｅｒ}について、実証実験で用いた物質（単軸性強誘電体結晶クロコン酸）で光子エネルギ１．５５ｅＶ（波長８００ｎｍ）において、α_{ｌａｓｅｒ}＝５２ｃｍ^－１＝０．００５２μｍ^－１である。ここで、深さの差が大きい層についてのテラヘルツ波８２の振幅を用いることで、深さの差の推定誤差が小さくなり、推定される吸収係数スペクトルの精度が高くなる。したがって、ここでは、図１２～図１５の例において、分析部１３０は、第７層及び第１層それぞれのＰ^－ _ｌ、及び、テラヘルツ波の振幅スペクトルから、レーザ光の吸収を考慮に入れたうえで解析を行う。

観測値とテラヘルツ帯の複素振幅との間には、以下の式１１で示す関係が成り立つ。

・・・（１１）

ここで、

・・・（１２）
は、第ｌ層で発生したテラヘルツ波形をフーリエ変換して得られたテラヘルツ帯の複素振幅スペクトルである。また、

・・・（１３）
は、テラヘルツ帯の複素屈折率スペクトルである。なお、便宜上、以後、式１３を「ｎ^～（ω）」と表記することがある。

ここで、テラヘルツ帯の複素屈折率スペクトルｎ^～（ω）について、以下の式１４が成り立つ。

・・・（１４）
なお、κ（ω）は、テラヘルツ帯の消衰係数スペクトルである。

式１１の左辺及び「ｚ_１－ｚ_７」は、因子分析の結果から求まる。また、α_{ｌａｓｅｒ}は、光学スペクトルの計測から独立に決定されているものとする。すると、式１１は、ｎ^～（ω）に関する方程式となり、ｎ^～（ω）が得られる。そして、式１４から、κ（ω）が求まる。そして、κ（ω）が求まれば、テラヘルツ帯の吸収係数スペクトルα（ω）は、一意に決定される。

なお、式１１の絶対値の方程式はα（ω）に関する方程式となり、以下の式１５で示すように解くことができる。したがって、テラヘルツ帯の吸収係数スペクトルα（ω）を一意に決定することができる。なお、式１５において、「ｌｏｇ」は自然対数である。なお、式１５において、「Ａｂｓ」は複素数の絶対値を得る演算子である。

・・・（１５）

このようにして、分析部１３０は、テラヘルツ帯の吸収係数スペクトルα（ω）を算出する。図１９は、図１２～図１５に示した例における、因子分析の結果を用いて推定されたテラヘルツ帯の吸収係数スペクトルを例示する図である。

［分析方法］
図５に示した分析方法について説明する。まず、分析システム１は、試料９０の照射位置ｋについて、ｋ＝１とする（ステップＳ１００）。具体的には、情報処理装置１００の分析装置制御部１１０は、照射位置ｋをｋ＝１に決定する。分析システム１は、照射位置ｋにフェムト秒レーザパルス２を照射させて、それによって発生するテラヘルツ波８０を観測する（ステップＳ１０２）。

具体的には、分析装置制御部１１０は、上述したように、分析装置１０（レーザ光源１２及び試料保持装置５０）を制御して、位置ｋに、フェムト秒レーザパルス２を照射させる。また、分析装置制御部１１０は、光学遅延ステージ４０を制御して、トリガパルス６が光伝導スイッチ３２に到達する時間を制御する。これにより、位置ｋで発生したテラヘルツ波８０の時間波形が得られる。そして、情報処理装置１００の取得部１２０は、観測されたテラヘルツ波８０を示す観測データを取得する（ステップＳ１０４）。このとき、取得部１２０は、観測データに照射位置ｋを対応付けて、観測データを記憶部１０４に記憶する。

分析システム１は、試料９０の照射位置ｋについて、ｋ＝Ｋであるか否かを判定する（ステップＳ１０６）。具体的には、情報処理装置１００の分析装置制御部１１０は、ｋ＝Ｋであるか否かを判定する。そして、ｋ＝Ｋでない場合（Ｓ１０６のＮＯ）、全ての照射位置ｋについて観測が終了していないので、分析装置制御部１１０は、ｋを１つインクリメントする（ステップＳ１０８）。そして、分析装置制御部１１０は、次の照射位置ｋについて、Ｓ１０２～Ｓ１０４の処理を繰り返す。これにより、レーザパルスがｘｙ面を走査されることとなる。そして、ｋ＝Ｋである場合（Ｓ１０６のＹＥＳ）、全ての照射位置ｋについて観測が終了する（ステップＳ１１２）。

分析部１３０は、得られた観測データに対して上述した因子分析を行う（ステップＳ１１４）。これにより、分析部１３０は、上述したように、分極行列Ｐ_ｋｌ及び波形ベクトルｖ_ｉを推定する（ステップＳ１１６）。また、分析部１３０は、上述したように、テラヘルツ帯の屈折率スペクトルを推定する（ステップＳ１１８）。また、分析部１３０は、上述したように、テラヘルツ帯の吸収係数スペクトルを推定する（ステップＳ１２０）。

さらに、分析結果出力部１４０は、分析結果を出力する（ステップＳ１２２）。具体的には、分析結果出力部１４０は、Ｓ１１６～Ｓ１２０で推定された結果を、インタフェース部１０８に出力させる。例えば、分析結果出力部１４０は、図１２～図１４に示す分析結果をインタフェース部１０８に表示させる。

［因子数の解析的決定方法］
ここで、因子数Ｆ（層数）をどのように設定するか、及びその信頼性の評価方法について述べる。解析的手法として、因子分析における「因子数」（Ｆ）を、パラメータとして設定する必要がある。因子数Ｆを増やせば増やすほどモデルが実験結果を再現するようになるが、あまりに因子数Ｆが大きくなりすぎると、ノイズを再現するだけになり、つまりアーティファクトを表現することとなる。したがって、因子分析においては、共分散行列の固有値をベイズ分析に基づいた赤池の情報基準（Akaike's Information Criterion：ＡＩＣ）などで分析し最適な因子数を決定する方法がある。

図２０は、実施の形態１にかかる、信頼できる因子数の決定方法を説明するための図である。本実施の形態にかかる方法では、共分散行列の固有値から適切な因子数Ｆを決定する。図２０は、因子数を十分大きくとって（図２０の例では因子数を２５個として）因子分析を行い、共分散行列の固有値を因子インデックスに対してプロットした結果を示す。この図２０に示したグラフは、スクリープロット（Scree Plot）と呼ばれる。図２０から、因子インデックス７以下と因子インデックス８以上では固有値の減少度合い（傾き）が異なることが分かる。因子分析においては、経験的に、この傾きの変化点が判断点になるとされている。これより、図２０の点線で示す直線に沿ったプロットに対応する因子インデックス８以上の因子分析の結果は、アーティファクトやノイズによるものと判断できる。このようにして、解析的に信頼できる層数（因子数）を設定することができる。これにより、試行錯誤による場合よりも、効率的にかつ適切に、因子数Ｆを決定することが可能となる。

なお、縦軸を対数スケールにしていることについては、テラヘルツ波形のデータ集合の持つ統計的性質による。
図２１は、テラヘルツ振幅の絶対値の３次元行列（ｘ，ｙ，ｔ依存性）のヒストグラムを示す図である。図２１から、データが指数分布をしていることが分かる。そのため、三次元像（テラヘルツ振幅のｘ，ｙ，ｔ依存性）の因子分析における固有値も指数分布をしたことがわかる。

［効果］
以下、本実施の形態にかかる技術の効果を、従来技術と比較しつつ説明する。特に、非特許文献１～４にかかる第二高調波発生（ＳＨＧ）による３次元ドメイン可視化手法との比較について述べる。

まず、ＳＨＧによるドメイン可視化の原理について述べる。強誘電体では、電気分極の向きにより二次の非線形光学係数（χ^（２））の符号が変化する。特に、単軸性の強誘電体（分極の向きが結晶の１つの軸方向に限られる強誘電体）の場合、この符号が変わる。すなわち、電気分極がχ^（２）と強く関連しているので、χ^（２）の変化を調べることで電気分極に関しての知見を得ることができる。実験的には、レーザ等の強いビームを強誘電体に照射し、第二高調波発生強度の場所依存性を測定することで、強誘電ドメインを可視化できる。その原理は、ビームが強誘電ドメイン内部に当たっている場合には、レーザ照射位置それぞれで発生したＳＨＧが同位相で重なるためにＳＨＧの強度が強くなり、一方、ビームが強誘電ドメイン境界に当たると、隣り合う強誘電ドメイン間でＳＨＧが逆位相であるためにＳＨＧが弱くなるということである。この原理は、基本波に対してＳＨＧが透過型であるか反射型であるかによらない。そして、強誘電体のＳＨＧ測定によって強誘電ドメインを可視化する研究は数多く行われている。

しかしながら、上記のようなＳＨＧの強度測定による強誘電ドメインイメージングでは、電気分極の大きさは分かるが、電気分極の符号を特定することは容易ではない。参照体からのＳＨＧと試料からのＳＨＧとを干渉させることで位相も含めて電気分極の符号を測定する手法があるが（非特許文献３）、装置や測定方法が複雑になる難点がある。また、試料の厚さをＳＨＧのコヒーレンス長より短くする必要があり、位相整合条件を満たすことが難しい試料の場合にはＳＨＧによる手法は有効ではない。また、ＳＨＧや差周波発生（ＤＦＧ）では、３次元可視化にある位相整合条件が満たされる必要がある（非特許文献１）。これに対し、本実施の形態にかかる方法では、時間領域で発生した電磁波（光）の振幅波形を観測できるため、このような制限条件は存在しない。

また、本実施の形態にかかる方法の利点として、長いコヒーレンス長による厚み不均一及び表面荒れに対するロバストさがあげられる。ＳＨＧを用いた三次元強誘電ドメイン可視化手法を実行するためには、コヒーレンス長を試料の厚さより長くすることが求められる。そのため、コヒーレンス長が試料厚さより長くなるようにするために、入射波長や偏光を制御するなど様々な工夫が行われている。しかしながら、この場合、使用することができるテンソル成分は制限され、結晶の対称性によっては位相整合条件を満たすことができないおそれがある。また、ＳＨＧは、試料の加工性が高い場合に有効である半面、波長の短さから散乱などの影響を受けやすい傾向がある。

これに対して、本実施の形態にかかる方法で用いるテラヘルツ電磁波の発生のコヒーレンス長は、ＳＨＧの場合のコヒーレンス長に比べて２桁程度長くなる傾向を持つ。例えば、ポンプ光とテラヘルツ電磁波との屈折率の差が１程度と大きいケースでも、透過配置の場合（入射方向と出射方向とが同じ方向の場合）のコヒーレンス長は、１ＴＨｚで１５０μｍとなる。このように、本実施の形態にかかる方法では、テラヘルツ電磁波の波長が長いため、厚みの不均一性や散乱の影響を受けにくい。さらに、上述したように、本実施の形態にかかる方法では、電場振幅を観測できる利点があり、したがって、非特許文献３にかかる技術で行われる干渉を用いずに電気分極の方向まで特定できる。

図２２は、本実施の形態にかかるドメインイメージング方法と他の光学的ドメインイメージング方法との比較を示す図である。図２２を用いて、テラヘルツ放射によるイメージング技術が他の光学的手法と比べてどのような特徴を有しているかについてさらに説明する。本実施の形態にかかる方法の１つの利点は、上述したように、分極方向の特定が容易であることである。

代表的な分光手法である蛍光や第二高調波発生においては、発生した光の「強度」をフォトダイオードなどで測定する。強誘電体での第二高調波発生（ＳＨＧ）の場合、ＳＨ光の強度は対称性の破れ度合いに比例する分極（Ｐ）の２乗に比例することが知られている。つまり、ＳＨＧでは、±Ｐのドメインを区別することは、ＳＨ光強度の測定のみでは可能でない。そこで、上述したように、参照物質を横に置き、その物質からのＳＨを試料からのＳＨと干渉させて計測することにより、±Ｐを区別する手法が開発されてきた。ただし、干渉信号取得のためには、試料端面をλ／１０レベルに光学研磨する必要がある。そのため、結晶の端面の光学研磨が容易な無機非線形光学結晶（ＬｉＮｂＯ_３など）にのみ、ＳＨＧの手法を適用可能であった。

これに対し、本実施の形態にかかる方法では、テラヘルツ放射の「振幅」は、対称性の破れ度（強誘電体の場合は分極Ｐ）に比例するので、光伝導スイッチや電気光学サンプリングによって符号を含めた振幅を測定することができる。そのため、干渉のための外部信号が不要である。また、発生するテラヘルツ波の波長は１００μｍ程度であるため、試料端面の厚みの不均一性に対してロバストである。例えば、１０μｍ程度の凸凹があってもテラヘルツ波の位相には１０％以下の変化率しか及ぼさない。

本実施の形態にかかるイメージング手法の特筆すべき技術優位性は、（１）非接触、（２）非破壊、（３）実デバイスへ使用可能（適用容易）、である。特に（３）については、外部電場による屈折率変調やピエゾ効果を用いる手法に比べて、強電場を印加しながら観測を行うことができるという長所がある。これに対し、第二高調波発生を用いる手法では、外部電場による屈折率変調が外部ＳＨＧ信号との干渉パターンに大きな影響を与えてしまうため、外部電場下での３次元ドメイン可視化は困難である。

また、本実施の形態にかかる方法では、放射されたテラヘルツ波形に対して統計的データ解析である因子分析を行う。これにより、３次元的に、分極値の推定を行うことができる。したがって、対象物質に関する光学スペクトルの情報の準備が最小限であっても、３次元ドメインの可視化を実現することが可能である。なお、上述したように、未知の対象物質の場合、光学スペクトル情報を予め取得することは困難であるおそれがある。したがって、本実施の形態にかかる分析システム１は、試料の結晶内部の３次元的な性質を効率的に可視化することが可能となる。なお、上述した非特許文献では、統計的データ解析を行って３次元ドメインの可視化を行うことについて、開示していない。したがって、３次元的なイメージング行うためには、対象物質に関する光学スペクトルの情報が必要である。したがって、試料の結晶内部の３次元的な性質を効率的に可視化することは困難である。

また、テラヘルツ波を試料に照射して、その波形に現れるエコー（反射波）から、試料の断面形状を解析する研究が行われている（例えば非特許文献６）。このような研究では、テラヘルツ波を試料に照射しているのであって、本実施の形態のように、レーザパルスを試料に照射してテラヘルツ波を試料から発生させているわけではない。また、このような研究では、試料の分極ドメインを可視化することは困難である。さらに、このような研究では、本実施の形態のように、統計的データ解析を行っているわけではない。したがって、３次元的なイメージング行うためには、対象物質に関する光学スペクトルの情報が必要である。したがって、試料の結晶内部の３次元的な性質を効率的に可視化することは困難である。

また、上述したように、本実施の形態にかかる分析システム１は、テラヘルツ帯の光学スペクトル（屈折率スペクトル及び吸収係数スペクトル）を抽出することも可能である。したがって、本実施の形態にかかる分析システム１は、未知の対象物質に対して効果的に分析を行うことが可能である。特に、テラヘルツ時間領域分光が困難な小さい物資への優位性がある。

また、本実施の形態にかかる分析システム１では、テラヘルツ波形測定によって絶対位相の情報が得られる。これにより、上述したような分析を行うことが可能となる。また、位相整合長の数倍（数ｍｍ）レベルのサイズの材料にも、本実施の形態にかかる方法は適用可能である。なお、強度測定では、位相の情報を得ることはできない。また、干渉法（非特許文献３）では、相対位相（０～３６０度）の情報のみ得られる。

また、本実施の形態にかかる方法は、厚い結晶に対するロバスト性も有している。ＳＨＧによる３次元イメージングでは、レーザを高開口数レンズで集光する必要がある。一方、本実施の形態にかかる方法では、低開口数のレンズで焦点深度を１ｍｍ程度にしても問題なく、深さ方向で全体に焦点が合った形で可視化を行うことが可能である。

（実施の形態２）
次に、実施の形態２について説明する。説明の明確化のため、以下の記載及び図面は、適宜、省略、及び簡略化がなされている。また、各図面において、同一の要素には同一の符号が付されており、必要に応じて重複説明は省略されている。実施の形態２では、多軸性強誘電体に対して本実施の形態にかかる方法を適用する例について説明する。なお、実施の形態２にかかるハードウェア構成及び機能構成は、図１及び図４に示したものと実質的に同様であるので、実施の形態２に特有の構成を除き、説明を省略する。

分極の向きが結晶の１つの軸方向に限られる強誘電体を単軸強誘電体という。例えば、実施の形態１では、分極の向きは、±方向の１つの軸方向である。一方、強誘電分極が結晶の複数の軸方向に向きうる強誘電体を多軸性強誘電体という。多軸性強誘電体においては、その非線形光学テンソルが変わるため、テラヘルツ波発生の偏光依存性が異なると考えられる。また、テラヘルツ波の発生機構が偏光により異なる場合や、レーザ又はテラヘルツ波に対する屈折率に偏光依存性があるために位相整合条件に偏光依存性が生じる場合は、テラヘルツ波形に入射偏光または検出偏光による変化が生じる。すなわち、深さ方向の同位置であっても、多軸性の異なる軸方向を向いたドメインについて発生したテラヘルツ波８０は異なる偏光依存性を示す。つまり、±方向以外の方向に対するドメインについて、異なる入射偏光依存性を示すテラヘルツ波８０が発生することになる。このように、多軸性の異なる軸方向を向いたドメインについて発生したテラヘルツ波８０に対してドメインの分布の推定及びベクトル分解を行うことによって、試料９０の３次元的な性質（ドメイン）をさらに精緻に分析することができる。

レーザパルス（ポンプパルス４）の偏光方向及び検出するテラヘルツ波８０の偏光方向を制御することで、多軸性のドメインにおいても、本実施の形態にかかる方法を適用できる。ここで、レーザパルス及びテラヘルツ波８０の偏光方向の設定値のインデックスをｑとし、その最大値をＱとする。つまり、レーザパルス及びテラヘルツ波８０の偏光方向の組み合わせは、Ｑ個である。また、Ｑは、偏光配置数ともいえる。なお、偏光配置ｑそれぞれがレーザパルス及びテラヘルツ波８０の偏光方向についてのどのような設定値に対応するかは、予め定められており、その対応付けは、情報処理装置１００によって予め記憶されている。

多軸性の強誘電体において生じ得るドメインの種類の数を、符号等の違いも考慮してＮ_{ｄｉｒｅｃｔｉｏｎｓ}とする。例えば、単軸性の場合は、ドメインの方向が＋方向及び－方向であるので、Ｎ_{ｄｉｒｅｃｔｉｏｎｓ}＝２である。また、例えば、平面内の二軸性強誘電体の場合は、Ｎ_{ｄｉｒｅｃｔｉｏｎｓ}＝２×２＝４となる。なお、ドメインの種類は、結晶の対称性を特徴づける要素であって、分極の方向だけでなく、例えば結晶構造の種類（α相又はβ相等）も含まれ得る。また、偏光方向の設定値ｑのそれぞれで２つのドメインの種類を区別できるので、偏光方向の設定値Ｑのテラヘルツ放射イメージングで区別できるドメインの種類数は２^Ｑである。したがって、Ｎ_{ｄｉｒｅｃｔｉｏｎｓ}種類のドメインを区別するためには、２^Ｑ≧Ｎ_{ｄｉｒｅｃｔｉｏｎｓ}とすればよい。または、冗長性を持たせるために、２^Ｑ＞Ｎ_{ｄｉｒｅｃｔｉｏｎｓ}を満たす偏光配置数Ｑでテラヘルツ波を観測してもよい。なお、本実施の形態では統計的データ解析を行っているため、冗長性を与えることで、より実際のドメイン構造を解析することができる。どのような偏光条件が最適かについては、結晶のもつ対称性から決定できる。

あるレーザスポット位置ｋにおいて、ある偏光配置ｑで観測した波形ベクトルを、以下の式１６のように表す。

・・・（１６）

上記の式１６を上述した因子分析に適用する場合、上記の式１６を、以下の式１７に示すような１次元ベクトルに変形する。以下の式１７で示す波形ベクトルｅ_ｋは、それぞれの偏光配置ｑで観測された波形ベクトルｅ_ｋｑを合成することによって得られた波形ベクトルである。なお、波形ベクトルｅ_ｋｑの成分がＭ個であるから、波形ベクトルｅ_ｋは、Ｍ×Ｑ個の成分で表され得る。

・・・（１７）

これを各レーザスポット位置ｋで測定し、実施の形態１（単軸強誘電体）の場合と同様の解析を行うことにより、全ての方向（ドメインの種類）を区別したドメイン可視化を実現することが可能である。なお、各因子に対応する分極ドメイン構造の深さ方向のインデックスについては、ぞれぞれの因子のテラヘルツ波形の位相スペクトルから決定できる。

図２３～図２５は、実施の形態２にかかる方法を説明するための図である。図２３～図２５を用いて、具体的なレーザパルス（ポンプパルス４）の偏光及び検出するテラヘルツ波８０の偏光方向の制御とその解析手順について述べる。テラヘルツ電磁波発生の偏光依存性の測定法としてよく用いられる方法は、入射光及び検出するテラヘルツ電磁波の偏光を固定し、試料を光軸周りで回転させる方法である。しかしながら、この方法は、試料が小さく、照射位置依存性も大きい場合には実験が困難であり、測定誤差も大きい。

そこで、実施の形態２では、ワイヤーグリッド偏光子２８Ａ，２８Ｂを用いたテラヘルツ波８０の偏光解析法と、１／２波長板２０による入射偏光の制御と、を組み合わせた方法を用いて、多軸性強誘電体の３次元的なドメイン可視化を実現する。この方法は、試料９０を回転又は移動させる必要がなく、試料９０の同じ位置にレーザを当てることができるため、照射位置のばらつきによる実験誤差を低減させる利点を有する。また、１／２波長板２０を自動回転させることで、測定の大部分を自動化できるメリットがある。

図２３は、実施の形態２にかかる、テラヘルツ波発生の偏光依存性測定の概念を説明するための図である。ポンプパルス４の光路上に配置された１／２波長板２０は、自動回転ステージ（図示せず）にセットされる。情報処理装置１００は、自動回転ステージを制御して、１／２波長板２０による入射光の偏光方向（偏光角θ_{ｌａｓｅｒ}）を制御する。１／２波長板２０を通過した直線偏光は、遅軸または速軸に対して折り返した直線偏光となる。

さらに、図２３の矢印Ａで示すように、実施の形態２にかかる分析システム１は、１／２波長板２０の回転角で制御した様々な偏光角θ_{ｌａｓｅｒ}の入射光について発生するテラヘルツ波８０の、水平方向（Ｘ軸方向）の偏光成分と垂直方向（Ｙ軸方向）の偏光成分とを、それぞれ測定する。そのために、テラヘルツ波８０の光路には、テラヘルツ放射源（試料９０）と光伝導スイッチ３２（テラヘルツ波８０の検出器）との間に、２つのワイヤーグリッド偏光子２８Ａ，２８Ｂが配置されている。ワイヤーグリッド偏光子２８は、ワイヤーに垂直な方向の偏光を透過し、ワイヤーに平行な方向の偏光を遮断する。

光伝導スイッチ３２の側のワイヤーグリッド偏光子２８Ｂ（ＷＧ２）は、光伝導アンテナ（光伝導スイッチ３２）の検出偏光を水平方向のみに制限する機能を有する。光伝導アンテナの間隙の方向は水平に固定されているとする。また、よく用いられるダイポール型の光伝導アンテナは、間隙の方向に平行なテラヘルツ波８０の電場成分を効率よく検出する。一方、間隙の方向に垂直な方向についての感度も、平行成分の感度の１００分の１程度あることが知られている。垂直成分の漏れは、このように成分測定上のノイズとなるため、垂直成分の漏れを完全にカットするために、ワイヤーグリッド偏光子２８Ｂ（ＷＧ２）を、検出器（光伝導スイッチ３２）と平行な水平方向を透過する配置で用いる。

テラヘルツ電磁波の偏光解析法について説明する。検出器としてダイポール型の光伝導スイッチ３２を用いる場合を考え、最も敏感に検出する方向をＸ軸方向とする。ただし、検出器に偏光依存性がない場合でも、ワイヤーグリッド偏光子２８ＢにおいてＸ軸方向の偏光のみを透過させるため、ベクトルマッピングを行うことが可能である。発生したテラヘルツ波８０の振幅ベクトルＥを、Ｅ＝（Ｅ_ｘ，Ｅ_ｙ）とする。テラヘルツ波８０の偏光をベクトルとして得るために、ワイヤーグリッド偏光子２８Ａ（ＷＧ１）の角度が±４５度の２つの配置それぞれでテラヘルツ波８０の波形の測定を行う。この場合、例えば、ワイヤーグリッド偏光子２８Ａは自動回転ステージ（図示せず）にセットされてもよい。そして、情報処理装置１００は、自動回転ステージを制御して、ワイヤーグリッド偏光子２８Ａの角度を制御する。なお、上述したように、ワイヤーグリッド偏光子２８Ｂ（ＷＧ２）の角度は０度に固定する。

また、ワイヤーグリッド偏光子２８Ａの角度が＋４５度，－４５度の配置で得られる信号を、それぞれ、Ｅ^{ｄｅｔｅｃｔｅｄ} _＋４５°，Ｅ^{ｄｅｔｅｃｔｅｄ} _－４５°とする。また、＋４５度，－４５度の配置のワイヤーグリッド偏光子２８Ａ（ＷＧ１）を透過した直後のテラヘルツ波の偏光ベクトル（振幅含む）を、それぞれ、Ｅ’_＋４５°，Ｅ’_－４５°とする。

また、＋４５度及び－４５度の方向の基底ベクトルを、それぞれ、以下の式１８及び式１９で示すように定める。

・・・（１８）

・・・（１９）

ここで、図２３の矢印Ｂで示すように、振幅ベクトルＥを＋４５度及び－４５度の方向に射影することを考える。このとき、ベクトルの射影の公式から、以下の式２０及び式２１が得られる。

・・・（２０）

・・・（２１）

ここで、ベクトルの射影の公式から、Ｘ軸方向及びＹ軸方向の基底ベクトルを、以下の式２２で表す。

・・・（２２）

このとき、ワイヤーグリッド偏光子２８Ｂ（ＷＧ２）の高透過率の偏光方向をｘ軸偏光としている場合、Ｅ^{ｄｅｔｅｃｔｅｄ} _＋４５°，Ｅ^{ｄｅｔｅｃｔｅｄ} _－４５°は、それぞれ以下の式２３，式２４のように表される。

・・・（２３）

・・・（２４）

式２３及び式２４で示される連立方程式を行列の形式で表すと、以下の式２５のように表される。

・・・（２５）

そして、両辺に係数行列の逆行列をかけると、以下の式２６で示される関係式が得られる。

・・・（２６）

上記の式２６より、発生したテラヘルツ波８０の振幅ベクトルＥの元々の水平偏光成分Ｅ_ｘ及び元々の垂直偏光成分Ｅ_ｙは、それぞれ、以下の式２７，式２８のように得られることが分かる。

・・・（２７）

・・・（２８）

したがって、ワイヤーグリッド偏光子２８Ａの角度が＋４５度，－４５度の配置で得られる信号Ｅ^{ｄｅｔｅｃｔｅｄ} _＋４５°，Ｅ^{ｄｅｔｅｃｔｅｄ} _－４５°から、テラヘルツ波８０の振幅ベクトル（電磁波ベクトル）Ｅ＝（Ｅ_ｘ，Ｅ_ｙ）を算出することができる。なお、上述した例では、ワイヤーグリッド偏光子２８Ａ（ＷＧ１）を±４５度に設定した場合について述べたが、ワイヤーグリッド偏光子２８Ａの角度は、±４５度に限定されるものではない。ベクトルの射影の公式を用いることで、任意のＷＧ１の角度の場合についても計算することができる。

図２４は、実施の形態２にかかる、各偏光配置で観測されたテラヘルツ波形を例示する図である。図２４は、１／２波長板２０の回転角で制御した偏光角θ_{ｌａｓｅｒ}を０度及び９０度としたときそれぞれについて、ワイヤーグリッド偏光子２８Ａを±４５度としたときに観測されるＥ^{ｄｅｔｅｃｔｅｄ} _＋４５°，Ｅ^{ｄｅｔｅｃｔｅｄ} _－４５°の波形を示している。

ここで、偏光角θ_{ｌａｓｅｒ}が９０度の場合にワイヤーグリッド偏光子２８Ａを±４５度としたときに観測される波形ベクトルＥ^{ｄｅｔｅｃｔｅｄ} _＋４５°，Ｅ^{ｄｅｔｅｃｔｅｄ} _－４５°を、それぞれＥ_１，Ｅ_２とする。また、偏光角θ_{ｌａｓｅｒ}が０度の場合にワイヤーグリッド偏光子２８Ａを±４５度としたときに観測される波形ベクトルＥ^{ｄｅｔｅｃｔｅｄ} _＋４５°，Ｅ^{ｄｅｔｅｃｔｅｄ} _－４５°を、それぞれＥ_３，Ｅ_４とする。また、Ｅ_１，Ｅ_２，Ｅ_３，Ｅ_４それぞれの偏光配置ｑのインデックスを、それぞれ、ｑ＝１，２，３，４とする（Ｑ＝４）。なお、図２４に示した偏光配置は、あくまでも例示である。偏光配置は、適宜設定され得る。

図２５は、実施の形態２にかかる因子分析の方法を説明するための図である。図２５に例示するように、実施の形態２では、観測された波形ベクトルＥ_１，Ｅ_２，Ｅ_３，Ｅ_４を、時間軸方向に並べる。つまり、Ｅ_１は式１７のｅ_ｋ１に対応し、Ｅ_２は式１７のｅ_ｋ２に対応し、Ｅ_３は式１７のｅ_ｋ３に対応し、Ｅ_４は式１７のｅ_ｋ４に対応する。そして、波形ベクトルＥ_１，Ｅ_２，Ｅ_３，Ｅ_４を時間軸方向に並べた波形ベクトルは、式１７のｅ_ｋに対応する。

実施の形態２では、分析部１３０は、各照射位置ｋについて、観測された波形ベクトルＥ_１，Ｅ_２，Ｅ_３，Ｅ_４を時間軸方向に並べて波形ベクトルｅ_ｋを生成する。そして、分析部１３０は、式６における波形ベクトルｖ_ｌの成分数をＭ×Ｑ個とした波形ベクトルｖ_ｌを共通因子として、実施の形態１と実質的に同様の因子分析を行う。なお、実施の形態２では、式６の行列Ｅ_{ＴＨｚ，ｋ}は、Ｋ×（Ｍ×Ｑ）の２次元行列となる。これにより、全ての方向の偏光に関するドメインを区別でき得る。したがって、実施の形態２にかかる分析システム１は、一軸性ドメイン構造に加えて、多軸性ドメイン構造を可視化することが可能となる。

図２６は、実施の形態２にかかる分析システム１によって実行される分析方法を示すフローチャートである。まず、図５のＳ１００と同様に、分析装置制御部１１０は、試料９０の照射位置ｋについて、ｋ＝１とする（ステップＳ２００）。さらに、分析装置制御部１１０は、偏光配置ｑについて、ｑ＝１とする（ステップＳ２０１）。

そして、分析システム１は、偏光配置ｑで照射位置ｋにフェムト秒レーザパルス２を照射させて、それによって発生するテラヘルツ波８０を観測する（ステップＳ２０２）。具体的には、分析装置制御部１１０は、上述したように、試料９０の位置ｋにフェムト秒レーザパルス２が照射されるように、試料保持装置５０を制御する。また、分析装置制御部１１０は、偏光配置ｑに対応する偏光方向の設定となるように、１／２波長板２０の回転角及びワイヤーグリッド偏光子２８Ａの回転角を制御する。そして、分析装置制御部１１０は、上述したように、レーザ光源１２を制御して、位置ｋに、フェムト秒レーザパルス２を照射させる。また、分析装置制御部１１０は、光学遅延ステージ４０を制御して、トリガパルス６が光伝導スイッチ３２に到達するタイミングを制御する。これにより、偏光配置ｑのときに位置ｋで発生したテラヘルツ波８０の時間波形が得られる。

ここで、例えば、各偏光配置ｑが１／２波長板２０の回転角で制御した偏光角θ_{ｌａｓｅｒ}とワイヤーグリッド偏光子２８Ａの角度との組み合わせによって設定される場合、情報処理装置１００は、各偏光配置ｑのときのこれらの組み合わせを記憶している。図２４の例では、情報処理装置１００は、偏光配置ｑ＝１のときの偏光角θ_{ｌａｓｅｒ}が９０度であること及びワイヤーグリッド偏光子２８Ａの角度が－４５度であることを記憶している。同様に、情報処理装置１００は、偏光配置ｑ＝２のときの偏光角θ_{ｌａｓｅｒ}が９０度であること及びワイヤーグリッド偏光子２８Ａの角度が＋４５度であることを記憶している。また、情報処理装置１００は、偏光配置ｑ＝３のときの偏光角θ_{ｌａｓｅｒ}が０度であること及びワイヤーグリッド偏光子２８Ａの角度が－４５度であることを記憶している。さらに、情報処理装置１００は、偏光配置ｑ＝４のときの偏光角θ_{ｌａｓｅｒ}が０度であること及びワイヤーグリッド偏光子２８Ａの角度が＋４５度であることを記憶している。したがって、例えば、ｑ＝１の場合、Ｓ２０２の処理において、分析装置制御部１１０は、偏光角θ_{ｌａｓｅｒ}が９０度となるように１／２波長板２０の回転角を制御する。また、分析装置制御部１１０は、ワイヤーグリッド偏光子２８Ａの角度が－４５度となるようにワイヤーグリッド偏光子２８Ａを制御する。

そして、情報処理装置１００の取得部１２０は、観測されたテラヘルツ波８０を示す観測データを取得する（ステップＳ２０３）。このとき、取得部１２０は、観測データに偏光配置ｑ及び照射位置ｋを対応付けて、観測データを記憶部１０４に記憶する。

分析システム１は、偏光配置ｑについて、ｑ＝Ｑであるか否かを判定する（ステップＳ２０４）。具体的には、情報処理装置１００の分析装置制御部１１０は、ｑ＝Ｑであるか否かを判定する。そして、ｑ＝Ｑでない場合（Ｓ２０４のＮＯ）、全ての偏光配置ｑについて観測が終了していないので、分析装置制御部１１０は、ｑを１つインクリメントする（ステップＳ２０５）。そして、分析装置制御部１１０は、次の偏光配置ｑについて、Ｓ２０２～Ｓ２０３の処理を繰り返す。そして、ｑ＝Ｑである場合（Ｓ２０４のＹＥＳ）、位置ｋについて、全ての偏光配置ｑでの観測が終了する。

分析システム１の分析装置制御部１１０は、試料９０の照射位置ｋについて、ｋ＝Ｋであるか否かを判定する（ステップＳ２０６）。そして、ｋ＝Ｋでない場合（Ｓ２０６のＮＯ）、全ての照射位置ｋについて観測が終了していないので、分析装置制御部１１０は、ｋを１つインクリメントする（ステップＳ２０８）。そして、分析装置制御部１１０は、次の照射位置ｋについて、Ｓ２０２～Ｓ２０５の処理を繰り返す。そして、ｋ＝Ｋである場合（Ｓ２０６のＹＥＳ）、全ての照射位置ｋについて、全ての偏光配置ｑでの観測が終了する（ステップＳ２１２）。

分析部１３０は、得られた観測データを用いて、各照射位置ｋについて、図２５に例示されるような波形を示す波形ベクトルｅ_ｋを生成する（ステップＳ２１４）。具体的には、分析部１３０は、照射位置ｋごとに、偏光配置ｑのインデックスの１からＱの順序で、観測データを時間軸方向に並べる。これにより、分析部１３０は、各照射位置ｋについての波形ベクトルｅ_ｋ（式１７）を生成する。

分析システム１は、得られた波形ベクトル（観測データ）を用いて、図５に示したＳ１１４～Ｓ１２２と実質的に同様の処理を行う（ステップＳ２２０）。具体的には、分析部１３０は、得られた波形ベクトルに対して、上述した因子分析を行う（Ｓ１１４）。これにより、分析部１３０は、分極行列Ｐ_ｋｌ及び波形ベクトルｖ_ｉを推定する（Ｓ１１６）。ここで、実施の形態２では、推定された分極行列Ｐ_ｋｌ及び波形ベクトルｖ_ｉには多軸性のドメインの特徴が示されている。すなわち、ｖ_ｉの符号およびｑ依存性から多軸性のドメインを区別することが可能である。ｖ_ｉに対応する分極行列Ｐ_ｋｌがそのドメインの分布を示している。

また、分析部１３０は、実施の形態１と同様にして、テラヘルツ帯の屈折率スペクトルを推定する（Ｓ１１８）。また、分析部１３０は、実施の形態１と同様にして、テラヘルツ帯の吸収係数スペクトルを推定する（Ｓ１２０）。さらに、分析結果出力部１４０は、分析結果を出力する（Ｓ１２２）。

実施の形態２にかかる分析システム１は、レーザパルス（レーザ光）の偏光方向及びテラヘルツ波の偏光方向の組み合わせである複数の偏光配置ｑのそれぞれについて、試料９０にレーザ光を入射させて、試料９０から出射したテラヘルツ波を検出する。そして、分析システム１は、照射位置ｋごとに、複数の偏光配置ｑそれぞれの観測データを時間軸方向に並べることで合成された波形ベクトル（ｅ_ｋ）を生成し、合成された波形ベクトルに対して統計的データ解析を行うことによって、ドメインの分布を推定する。これにより、多軸性の強誘電体において生じ得る複数のドメインの種類について、ドメインの分布を推定する。したがって、複数のドメインの方向や種類を持つ物質に対して、複数のドメインを区別したドメインの可視化を行うことができる。

（実施の形態３）
次に、実施の形態３について説明する。説明の明確化のため、以下の記載及び図面は、適宜、省略、及び簡略化がなされている。また、各図面において、同一の要素には同一の符号が付されており、必要に応じて重複説明は省略されている。

図２７は、実施の形態３にかかる光学配置の例を示す図である。上述した実施の形態では、図２７の矢印Ａで示すように、フェムト秒レーザパルス２が試料９０に入射する入射方向と試料９０からテラヘルツ波８０が出射する出射方向とが同一方向である（透過配置）。しかしながら、フェムト秒レーザパルス２の入射方向とテラヘルツ波８０の出射方向とは、同一方向である必要はない。

例えば、図２７の矢印Ｂ～Ｄで示すように、フェムト秒レーザパルス２の入射方向とテラヘルツ波８０の出射方向とは、同一方向でなくてもよい。図２７の矢印Ｂは、フェムト秒レーザパルス２が試料９０に垂直に入射して、テラヘルツ波８０が垂直に反射する配置（反射配置）を示す。また、図２７の矢印Ｃは、フェムト秒レーザパルス２が試料９０に垂直に入射して、テラヘルツ波８０が斜め方向に反射する配置（斜出射検出配置）を示す。また、図２７の矢印Ｄは、フェムト秒レーザパルス２が試料９０に斜め方向から入射して、テラヘルツ波８０が斜め方向に反射する配置（斜出射検出配置）を示す。このような場合であっても、上述した実施の形態と同様に、３次元的なドメインの可視化が可能である。

なお、このように、フェムト秒レーザパルス２の入射方向とテラヘルツ波８０の出射方向とが同一方向でない場合、深さ方向の各層（深さ位置ｌそれぞれ）で発生するテラヘルツ波８２の時間遅延がより大きくなるため、深さ方向の分解能が向上する可能性がある。テラヘルツ波の最大周波数ｆ_ｍａｘ＝３ＴＨｚとすると、その波長は１００μｍである。ここで、テラヘルツ波観測ではπ／６程度の位相差を検知できることから、深さ方向の分解能としては、１００／１２μｍ～１０μｍを実現できる（但し「～」は近似記号）。反射配置は、さらに、テラヘルツ波又はレーザ光に対して透過率が低い場合に有効である。

（実施の形態４）
次に、実施の形態４について説明する。説明の明確化のため、以下の記載及び図面は、適宜、省略、及び簡略化がなされている。また、各図面において、同一の要素には同一の符号が付されており、必要に応じて重複説明は省略されている。

上述した実施の形態では、Δｎ（ω）＞０（つまり、式８よりｎ（ω）＞ｎ_ｇ）であるとしていたが、試料によっては、Δｎ（ω）＜０（つまり、ｎ（ω）＜ｎ_ｇ）である場合もある。本発明は、Δｎ（ω）＜０の試料についても適用可能である。以下、Δｎ（ω）＜０の場合も考慮した方法について述べる。具体的には、実施の形態４では、レーザ波長での群屈折率で決まる、レーザ光の包括線の伝搬速度（ｃ／ｎ_ｇ）とテラヘルツ波の伝搬速度（ｃ／ｎ（ω））の大小によって本手法がどのような影響を受けるかについて述べる。

上述した実施の形態にかかる実証に用いた試料（クロコン酸結晶）においては、観測されるテラヘルツ波の帯域全体でｎ（ω）＞ｎ_ｇの条件が満たされていることは、テラヘルツ時間領域分光の結果からすでに明らかになっていた。しかし、多くの物質で屈折率の分散の傾向はｎ（ω）＞ｎ_ｇであるものの、この関係式が一般の未知の物質において成り立っているとは限らない。

ここで、上述したように、ｎ（ω）＞ｎ_ｇである物質の場合、位相スペクトルθ（ω）の曲線の傾きがマイナス側に大きいテラヘルツ波の発生した深さ位置が、レーザ入射側の深さ位置に対応する。一方、ｎ（ω）＜ｎ_ｇである物質の場合は、位相スペクトルθ（ω）の曲線の傾きがプラス側に大きいテラヘルツ波の発生した深さ位置が、レーザ入射側の深さ位置に対応する。すなわち、時間軸上では、レーザ入射側の面（深さ位置）で発生したテラヘルツ波の方が、レーザパルスよりもより早く伝搬し、レーザパルスの信号よりもマイナス時刻側にテラヘルツ波の信号が出現する。このようにして、ｎ（ω）＜ｎ_ｇの試料の場合でも、位相曲線の傾きから、因子成分ｂと深さ位置ｌとの対応を判定することができる。

図２８は、実施の形態４にかかる方法を説明するための図である。図２８を用いて、ｎ（ω）とｎ_ｇの大小関係を概算評価する方法を述べる。ここでは、レーザ光に対して試料は十分透明であるとする。まず、レーザパルス（フェムト秒レーザパルス２）に対する群屈折率ｎ_ｇは、レーザ波長域での屈折率スペクトルから計算することができ、これは本発明での解析では既知としてよい。一方、テラヘルツ帯の屈折率スペクトルｎ（ω）は未知である。

そこで、可視・テラヘルツ帯の屈折率情報が知られており、観測されるテラヘルツ周波数帯で位相整合条件がよく満たされている結晶を参照物質として、同条件でテラヘルツ波形を測定することを考える。ここでは、チタンサファイアモードロックレーザの波長（８００ｎｍ）で非線形光学効果によるテラヘルツ波発生の位相整合条件がよく満たされている結晶を参照物質として用いる。チタンサファイアモードロックレーザの波長（８００ｎｍ）では、ＺｎＴｅ結晶がよく使用される。参照物質はｎ_ｇ ^ｒｅｆ～ｎ_ＴＨｚ ^ｒｅｆで位相整合条件がよく満たされているとする（但し「～」は近似記号）。ここで、ｎ_ＴＨｚ ^ｒｅｆは、テラヘルツ放射の参照物質のテラヘルツ帯の屈折率の平均値である。この場合、観測されるテラヘルツ波はほぼ全周波数で同位相となるため、シャープなピークをとることが知られている。

図２８は、レーザパルス及びテラヘルツ波が伝播した際のある時刻での位置を示している。図２８において、横軸は光軸方向を示す。図２８の（ａ）は、物体（試料）がない場合（つまり真空の場合）の、フェムト秒レーザパルス２の伝播を示す。図２８の（ａ）は、ある時刻におけるフェムト秒レーザパルス２の位置を示す。また、図２８の（ｂ）は、フェムト秒レーザパルス２、及び、位相整合条件が良い参照物質で発生したテラヘルツ波の伝播を示す。図２８の（ｃ）は、フェムト秒レーザパルス２、及び、ｎ（ω）＞ｎ_ｇの試料で発生したテラヘルツ波の伝播を示す。図２８の（ｄ）は、フェムト秒レーザパルス２、及び、ｎ（ω）＜ｎ_ｇの試料で発生したテラヘルツ波の伝播を示す。

ｎ（ω）とｎ_ｇの大小の評価方法について述べる。参照物質の厚さをｄ^ｒｅｆとすると、真空を伝搬する場合と比べ、光路長は、以下に示す式２９のように増加している。ここで、ｎ_ｇ ^ｒｅｆは、テラヘルツ放射の参照物質におけるレーザ波長での群屈折率である。この群屈折率ｎ_ｇ ^ｒｅｆは既知であるとする。

・・・（２９）
したがって、以下の式３０で示す時間遅延が生じていると看做される。δｔ^ｒｅｆは、真空を伝搬した場合と比べた、参照物質伝搬時の時間遅延を示す。

・・・（３０）
これにより、後述する、試料伝搬時のレーザ光及びテラヘルツ波の時間遅延の、基準点が定まる。つまり、図２８において、ｃδｔ^ｒｅｆから、矢印Ａで示す遅延の基準点Ｐが定まる。

一方、試料のレーザ入射側で発生したテラヘルツ波が試料内を伝搬し、放射される場合の時間遅延は、以下の式３１で表される。式３１は、真空を伝搬した場合（基準点Ｐに対応）と比べた、試料伝搬時のテラヘルツ波の時間遅延を示す。式３１において、ｎは観測テラヘルツ帯域での屈折率の平均値である。

・・・（３１）

また、テラヘルツ波出射側でテラヘルツ波が発生した場合の時間遅延は、以下の式３２で表される。式３２は、真空を伝搬した場合（基準点Ｐに対応）と比べた、試料伝搬時のレーザ光の時間遅延を示す。

・・・（３２）

図２８の（ｃ）で示すように、実験で観測されたテラヘルツ波の第一のメインピーク（最も早く観測点に到達したピーク）の時刻がδｔ_{ｌａｓｅｒ} ^ｓａｍから予想される位置にある場合、レーザ光がテラヘルツ波よりも早く伝搬している。この場合、この試料では、ｎ（ω）＞ｎ_ｇとみなせる。一方、図２８の（ｄ）で示すように、実験で観測されたテラヘルツ波の第一のメインピークの時刻がδｔ_{ｌａｓｅｒ} ^ｓａｍからの予測位置よりも早い場合、この関係式ｎ（ω）＞ｎ_ｇは成り立っていない。したがって、この試料では、ｎ（ω）＜ｎ_ｇと予測される。なお、図２８の（ｃ）及び（ｄ）の黒丸は、試料においてテラヘルツ波の第一のメインピークが発生したと予想される点を示す。つまり、ｎ（ω）＞ｎ_ｇの試料では、試料のテラヘルツ波出射側で発生したテラヘルツ波が、最も早く観測点に到達すると予想される。一方、ｎ（ω）＜ｎ_ｇの試料では、試料のレーザ入射側で発生したテラヘルツ波が、最も早く観測点に到達すると予想される。このようにして、ｎ（ω）とｎ_ｇの大小関係は、試料からのテラヘルツ波形の第一のメインピークの位置を、参照物質からの波形のメインピークの位置と比較することにより、判定することができる。

（変形例）
なお、本発明は上記実施の形態に限られたものではなく、趣旨を逸脱しない範囲で適宜変更することが可能である。例えば、上述した実施の形態では、試料が強誘電体であるとしたが、このような構成に限られない。本実施の形態は、強磁性体、焦電体又は非線形光学結晶にも適用可能である。また、上述した実施の形態では、試料の電気分極の３次元ドメインを可視化するとしたが、このような構成に限られない。例えば、試料が強磁性体である場合、磁気分極（磁化）の３次元ドメインを可視化するようにしてもよい。

また、上述した実施の形態では、試料の３次元的な分極の分布及び試料の各深さ位置で発生したテラヘルツ波の波形を推定するために、因子分析を用いているが、このような構成に限定されない。例えば、動的因子分析（Dynamic factor analysis）を用いてもよい。すなわち、試料の厚さが一様でない場合は、それぞれの照射位置で観測されるテラヘルツ波に時間シフトが生じる。この場合、動的因子分析（Dynamic factor analysis）を適用することで、厚さの揺らぎによる時間軸上でのシフトを考慮して、上記の推定を行うことができる。

また、上述した実施の形態では、試料９０にフェムト秒レーザパルス２等のレーザパルスを照射することによってテラヘルツ波８０を発生させるとしたが、このような構成に限られない。テラヘルツ波８０を発生させることが可能な任意のレーザ光を照射してもよい。

また、上述したフローチャートにおいて、各処理（ステップ）の順序は、適宜、変更可能である。また、複数ある処理（ステップ）のうちの１つ以上は、省略されてもよい。例えば、図５のＳ１１８，Ｓ１２０の処理の少なくとも一方はなくてもよい。また、図２６においては、ある照射位置ｋについて全ての偏光配置ｑ（ｑ＝１～Ｑ）で観測を行ってから、次の照射位置ｋ＋１で同様の観測を行うようにしているが、このような順序に限られない。ある偏光配置ｑについて全ての照射位置ｋで観測を行ってから、次の偏光配置ｑ＋１で同様の観測を行うようにしてもよい。つまり、図２６において、ｋとｑとを入れ替えてもよい。

１ 分析システム
２ フェムト秒レーザパルス
４ ポンプパルス
６ トリガパルス
１０ 分析装置
１２ レーザ光源
１４ ビームスプリッタ
１６ オプティカルチョッパ
１８ 反射型可変ＮＤフィルタ
２０ １／２波長板
２２ 近赤外対物レンズ
２４，３０ 軸外し放物面鏡
２６ 黒色低密度ポリエチレンフィルム
２８ ワイヤーグリッド偏光子
３２ 光伝導スイッチ
４０ 光学遅延ステージ
４２ 可動鏡
４４ レンズ
４６ 電流電圧変換器
４８ ロックイン検出器
５０ 試料保持装置
５２ 試料保持部
５４，５６ 可動ステージ
８０，８２ テラヘルツ波
９０ 試料
９０ａ 入射面
９０ｂ 出射面
１００ 情報処理装置
１１０ 分析装置制御部
１２０ 取得部
１３０ 分析部
１４０ 分析結果出力部

Claims (16)

1. 試料の第１の面に沿ってレーザ光を走査して前記レーザ光を前記第１の面の側から前記試料に入射させ、前記試料から出射したテラヘルツ波を検出する分析装置と、
   前記テラヘルツ波の観測された振幅の時間的な推移を示す観測データを取得する取得部と、
   取得された前記観測データに対して統計的データ解析を行うことによって、前記試料の、前記レーザ光の入射方向に沿った方向である深さ方向の複数の深さ位置それぞれにおける、前記レーザ光を走査した位置に対応する照射位置に関するドメインの分布を推定する分析部と、
   を有する分析システム。
2. 前記取得部は、前記試料の前記複数の深さ位置それぞれで発生し前記試料内を伝搬した後に前記試料の外に出射された前記テラヘルツ波に関する前記観測データを取得し、
   前記分析部は、取得された前記観測データに対して統計的データ解析を行うことによって、前記取得部で取得された前記観測データに対応する前記テラヘルツ波の波形を推定する、
   請求項１に記載の分析システム。
3. 前記分析部は、前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する位相スペクトルを用いて、推定された前記ドメインの分布に対応する前記深さ位置を推定する、
   請求項２に記載の分析システム。
4. 前記分析部は、取得された前記観測データに対して因子分析を行うことによって、前記ドメインの分布を推定する、
   請求項１から３のいずれか１項に記載の分析システム。
5. 前記分析部は、前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する位相スペクトルを用いて、テラヘルツ帯の屈折率スペクトルを推定する、
   請求項２から４のいずれか１項に記載の分析システム。
6. 前記分析部は、前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する複素振幅スペクトルを用いて、テラヘルツ帯の吸収係数スペクトルを推定する、
   請求項２から５のいずれか１項に記載の分析システム。
7. 前記分析装置は、前記レーザ光の偏光方向及び前記テラヘルツ波の偏光方向の組み合わせである複数の偏光配置のそれぞれについて、前記試料にレーザ光を入射させて、前記試料から出射した前記テラヘルツ波を検出し、
   前記分析部は、前記照射位置ごとに、複数の前記偏光配置それぞれの前記観測データを時間軸方向に並べることで合成された波形ベクトルを生成し、合成された波形ベクトルに対して統計的データ解析を行うことによって、前記ドメインの分布を推定する、
   請求項１から６のいずれか１項に記載の分析システム。
8. 試料の第１の面に沿ってレーザ光を走査して前記レーザ光を前記第１の面の側から前記試料に入射させ、前記試料から出射したテラヘルツ波を検出する分析装置を制御する分析装置制御部と、
   前記テラヘルツ波の観測された振幅の時間的な推移を示す観測データを取得する取得部と、
   取得された前記観測データに対して統計的データ解析を行うことによって、前記試料の、前記レーザ光の入射方向に沿った方向である深さ方向の複数の深さ位置それぞれにおける、前記レーザ光を走査した位置に対応する照射位置に関するドメインの分布を推定する分析部と、
   を有する情報処理装置。
9. 試料の第１の面に沿ってレーザ光を走査して前記レーザ光を前記第１の面の側から前記試料に入射させ、前記試料から出射したテラヘルツ波を検出し、
   前記テラヘルツ波の観測された振幅の時間的な推移を示す観測データを取得し、
   取得された前記観測データに対して統計的データ解析を行うことによって、前記試料の、前記レーザ光の入射方向に沿った方向である深さ方向の複数の深さ位置それぞれにおける、前記レーザ光を走査した位置に対応する照射位置に関するドメインの分布を推定する、
   分析方法。
10. 前記試料の前記複数の深さ位置それぞれで発生し前記試料内を伝搬した後に前記試料の外に出射された前記テラヘルツ波に関する前記観測データを取得し、
    取得された前記観測データに対して統計的データ解析を行うことによって、取得された前記観測データに対応する前記テラヘルツ波の波形を推定する、
    請求項９に記載の分析方法。
11. 前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する位相スペクトルを用いて、推定された前記ドメインの分布に対応する前記深さ位置を推定する、
    請求項１０に記載の分析方法。
12. 取得された前記観測データに対して因子分析を行うことによって、前記ドメインの分布を推定する、
    請求項９から１１のいずれか１項に記載の分析方法。
13. 前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する位相スペクトルを用いて、テラヘルツ帯の屈折率スペクトルを推定する、
    請求項１０から１２のいずれか１項に記載の分析方法。
14. 前記複数の深さ位置それぞれで発生したテラヘルツ波の波形に対応する複素振幅スペクトルを用いて、テラヘルツ帯の吸収係数スペクトルを推定する、
    請求項１０から１３のいずれか１項に記載の分析方法。
15. 前記レーザ光の偏光方向及び前記テラヘルツ波の偏光方向の組み合わせである複数の偏光配置のそれぞれについて、前記試料にレーザ光を入射させて、前記試料から出射した前記テラヘルツ波を検出し、
    前記照射位置ごとに、複数の前記偏光配置それぞれの前記観測データを時間軸方向に並べることで合成された波形ベクトルを生成し、合成された波形ベクトルに対して統計的データ解析を行うことによって、前記ドメインの分布を推定する、
    請求項９から１４のいずれか１項に記載の分析方法。
16. 試料の第１の面に沿ってレーザ光を走査して前記レーザ光を前記第１の面の側から前記試料に入射させ、前記試料から出射したテラヘルツ波を検出する分析装置を制御するステップと、
    前記テラヘルツ波の観測された振幅の時間的な推移を示す観測データを取得するステップと、
    取得された前記観測データに対して統計的データ解析を行うことによって、前記試料の、前記レーザ光の入射方向に沿った方向である深さ方向の複数の深さ位置それぞれにおける、前記レーザ光を走査した位置に対応する照射位置に関するドメインの分布を推定するステップと、
    をコンピュータに実行させるプログラム。

Images