特許 7571160 サンプル特性評価のための超高分解能時間領域分光法

(19)【発行国】日本国特許庁(JP)

(12)【公報種別】特許公報(B2)

(11)【特許番号】

(24)【登録日】2024-10-11

(45)【発行日】2024-10-22

(54)【発明の名称】サンプル特性評価のための超高分解能時間領域分光法

(51)【国際特許分類】

   G01N 21/3586 20140101AFI20241015BHJP

【ＦＩ】

G01N21/3586

【請求項の数】 11

(21)【出願番号】P 2022573798

(86)(22)【出願日】2021-02-11

(65)【公表番号】

(43)【公表日】2023-03-30

(86)【国際出願番号】 EP2021053298

(87)【国際公開番号】W WO2021160725

(87)【国際公開日】2021-08-19

【審査請求日】2023-11-27

(31)【優先権主張番号】20305128.9

(32)【優先日】2020-02-11

(33)【優先権主張国・地域又は機関】EP

(73)【特許権者】

【識別番号】311016455

【氏名又は名称】サントル ナシオナル ドゥ ラ ルシェルシェ シアンティフィク

【氏名又は名称原語表記】ＣＥＮＴＲＥ ＮＡＴＩＯＮＡＬ ＤＥ ＬＡ ＲＥＣＨＥＲＣＨＥ ＳＣＩＥＮＴＩＦＩＱＵＥ

(73)【特許権者】

【識別番号】518202378

【氏名又は名称】サントラル・リール・アンスティテュ

(73)【特許権者】

【識別番号】522320006

【氏名又は名称】ウニベルシテ・ドゥ・リール

(73)【特許権者】

【識別番号】522222489

【氏名又は名称】ウニベルシテ・ポリテクニック・オー－ドゥ－フランス

(73)【特許権者】

【識別番号】522320017

【氏名又は名称】ジュニア

(74)【代理人】

【識別番号】110001173

【氏名又は名称】弁理士法人川口國際特許事務所

(72)【発明者】

【氏名】ペレッティ，ロマン

【審査官】井上 徹

(56)【参考文献】

【文献】国際公開第２０１４／０３４０８５（ＷＯ，Ａ１）

【文献】米国特許出願公開第２０１７／０３０７５２０（ＵＳ，Ａ１）

(58)【調査した分野】(Int.Cl.，ＤＢ名)

Ｇ０１Ｎ ２１／００－２１／０１

Ｇ０１Ｎ ２１／１７－２１／６１

Ｇ０１Ｊ ３／００－ ４／０４

Ｇ０１Ｊ ７／００－ ９／０４

ＪＳＴＰｌｕｓ／ＪＭＥＤＰｌｕｓ／ＪＳＴ７５８０（ＪＤｒｅａｍＩＩＩ）

(57)【特許請求の範囲】

【請求項1】

サンプルの物理パラメータ（ｐｉ）の組を特定する方法であって、
－ Ａ 測定サンプル時間トレースＥｓ（ｔ）を取得するステップであって、前記測定サンプル時間トレースＥｓ（ｔ）は、時間領域分光法により、期間Ｔを有し、且つコム周波数を示す電磁パルスを周期的に発出する励起ビーム（ＥＢ）によってサンプル（Ｓ）を照明し、及び前記サンプルに由来する電磁場をコヒーレント検出によって時間の関数として検出することによって取得されており、サンプル時間トレースが測定される持続時間は、ｔｍａｘであり、ここで、ｔｍａｘ＜Ｔである、ステップ、
－ Ｂ 測定参照時間トレースＥｒｅｆ（ｔ）を取得するステップであって、前記測定参照時間トレースＥｒｅｆ（ｔ）は、ステップＡと同じ条件において、ただし前記サンプルの存在なしでの照明及び検出によって取得されている、ステップ、
－ Ｃ 前記期間Ｔにおいて広がり、且つ測定が行われなかった時点にゼロ値を付与することによって得られる、Ｅｒｅｆ０（ｔ）と呼ばれる拡張参照時間トレースを特定し、及びＴに等しい時間ウィンドウについて計算される、前記拡張参照時間トレースの離散フーリエ変換

【数1】

を特定するステップ、
－ Ｄ 前記拡張参照時間トレースの前記フーリエ変換

【数2】

及び前記サンプルの物理挙動モデルから、前記物理パラメータ（ｐｉ）の組に依存する、サンプル周波数モデル

【数3】

と呼ばれる、周波数領域での前記サンプルのインパルス応答のモデリングを特定するステップ、
－ Ｅ 前記物理パラメータ（ｐｉ）の組に最適化アルゴリズムを適用するステップであって、
－ Ｅ１ 物理パラメータ（ｐｉ）を初期化するサブステップ、
－ 以下のサブステップ：
－ Ｅ２ 推定サンプル時間トレースＥ_ｅｓｔ｛ｐ_ｉ｝（ｔ）と呼ばれる、前記サンプル周波数モデル

【数4】

の逆離散フーリエ変換を計算するサブステップ、
－ Ｅ３ 関数ｆに従って、前記測定サンプル時間トレースＥｓ（ｔ）と推定時間トレースＥｅｓｔ（ｔ）との間の差からエラー関数（ε_ｅｒ｛ｐｉ｝）を計算するサブステップであって、ｆは、形態学的距離を定義する関数であり、

【数5】

である、サブステップ
を、前記エラー関数を最小化する物理パラメータの値（ｐｉ_ｏｐｔ）の組を得るまで繰り返し実行するサブステップ
を含むステップ
を含む方法。

【請求項2】

前記励起ビーム（ＥＢ）は、ＴＨｚ領域内にあり、１００ＧＨｚ～３０ＴＨｚに含まれる周波数を有する、請求項１に記載の方法。

【請求項3】

前記サンプル時間トレースが測定される持続時間ｔｍａｘは、前記測定参照時間トレースのエネルギーの９５％超を含むように選択される、請求項１又は２に記載の方法。

【請求項4】

前記サンプル周波数モデル

【数6】

は、フーリエ変換Ｅｒｅｆ０（ω）に、サンプル挙動を特徴付ける伝達関数Ｔ（ω）を乗じることを含む、請求項１～３の何れか一項に記載の方法。

【請求項5】

前記伝達関数Ｔ（ω）は、複素屈折率ｎ（ω）に依存する、請求項４に記載の方法。

【請求項6】

誘電率ε（ω）と呼ばれる前記複素屈折率の二乗は、各スペクトル線に関するドルーデ・ローレンツモデルに追従し、スペクトル線は、３つのパラメータ、振幅（Ｍ）、減衰率（γ）と呼ばれる幅及び中心周波数（ω０）の組によって特徴付けられる、請求項５に記載の方法。

【請求項7】

前記エラー関数は、

【数7】

として定義される、請求項１～６の何れか一項に記載の方法。

【請求項8】

サンプル（Ｓ）を特性評価するための特性評価装置（５０）であって、
－ 測定サンプル時間トレースＥｓ（ｔ）及び測定参照時間トレースＥｒｅｆ（ｔ）を記憶するメモリ（ＭＥＭ）であって、前記測定サンプル時間トレースＥｓ（ｔ）は、時間領域分光法により、期間Ｔを有し、且つコム周波数を示す電磁パルスを周期的に発出する励起ビーム（ＥＢ）によって前記サンプル（Ｓ）を照明し、及び前記サンプルに由来する電磁場をコヒーレント検出によって時間の関数として検出することによって取得されており、サンプル時間トレースが測定される持続時間は、ｔｍａｘであり、ここで、ｔｍａｘ＜Ｔであり、
前記測定参照時間トレースＥｒｅｆ（ｔ）は、測定サンプル時間トレースＥｓ（ｔ）と同じ条件において、ただし前記サンプルの存在なしでの照明及び検出よって取得されている、メモリ（ＭＥＭ）、
－ 処理ユニット（ＰＵ）であって、
－ 前記期間Ｔにおいて広がり、且つ測定が行われなかった時点にゼロ値を付与することによって得られる、Ｅｒｅｆ０（ｔ）と呼ばれる拡張参照時間トレースを特定し、及びＴに等しい時間ウィンドウについて計算される、前記拡張参照時間トレースの離散フーリエ変換

【数8】

を特定すること、
－ 前記拡張参照時間トレースの前記フーリエ変換Ｅｒｅｆ０（ω）及び前記サンプルの物理挙動モデルから、物理パラメータ（ｐｉ）の組に依存する、サンプル周波数モデル

【数9】

と呼ばれる、周波数領域での前記サンプルのインパルス応答のモデリングを特定すること、
－ 前記物理パラメータ（ｐｉ）に最適化アルゴリズムを適用することであって、
－ 物理パラメータを初期化するステップ、
－ 以下のサブステップ：
－ 推定サンプル時間トレースＥ_ｅｓｔ｛ｐ_ｉ｝（ｔ）と呼ばれる、前記サンプル周波数モデル

【数10】

の逆離散フーリエ変換を計算するサブステップ、
－ 関数ｆに従って、前記測定サンプル時間トレースＥｓ（ｔ）と推定時間トレースＥｅｓｔ（ｔ）との間の差からエラー関数（ε_ｅｒ｛ｐｉ｝）を計算するサブステップであって、ｆは、形態学的距離を定義する関数であり、

【数11】

である、サブステップ
を、前記エラー関数を最小化する物理パラメータの値（ｐ_ｏｐｔｉ）の組を得るまで繰り返し実行するステップ
を含む、適用すること
を行うように構成された処理ユニット（ＰＵ）
を含む特性評価装置（５０）。

【請求項9】

分光光度計（Ｓｐｅｃｔｒｏ）であって、
－ 請求項８に記載の特性評価装置（５０）、
－ 測定装置（ＭｅＤ）であって、
－ サンプル（Ｓ）を前記励起ビーム（ＥＢ）によって照明するように構成された光源（ＬＳ）、
－ 前記測定サンプル時間トレースＥｓ（ｔ）及び前記測定参照時間トレースＥｒｅｆ（ｔ）を検出するように構成された検出器（Ｄ）
を含む測定装置（ＭｅＤ）
を含む分光光度計（Ｓｐｅｃｔｒｏ）。

【請求項10】

請求項１～７の何れか一項に記載のステップを実施するように適合されたコンピュータプログラム。

【請求項11】

請求項１０に記載のコンピュータプログラムを組み込んだコンピュータ可読媒体。

【発明の詳細な説明】

【技術分野】

【0001】

本発明は、概して、時間領域分光法（ＴＤＳ）により、サンプルの物理パラメータを特定する方法に関する。特に、本発明は、サンプルの測定時間トレースをフィッティングし、且つそのようにして物理パラメータを抽出する方法に関する。

【背景技術】

【0002】

時間領域分光法は、今日、テラヘルツ時間領域（ＴＨｚ－ＴＤＳ）において、すなわち０．１ＴＨｚ～３０ＴＨｚの周波数の励起ビームを使用することによって実行される定着した技術である。ＴＨｚ－ＴＤＳは、半導体、強誘電体、超導電体、液体、気体、生体分子、炭水化物等の分子結晶、バンドパスフィルタ、マイクロ流体回路構成に埋め込まれたメタサーフェス等、様々な材料を調査できることが証明されている。

【0003】

測定関数が干渉法を通した時間領域データの自己相関であるフーリエ変換赤外分光計とは対照的に、ＴＨｚ－ＴＤＳは、ＴＨｚ周波数範囲内の電界を直接測定することに基づく。

【0004】

一般的なＴＨｚ－ＴＤＳ装置の動作原理を図１に示す。ＴＨｚ励起ビームＥＢは、光伝導性ＴＨｚアンテナ又は非線形結晶ＮＬにより、フェムト秒レーザＬによって生成される近赤外パルスの光整流効果によって発出される（図１に示す）。励起ビームＥＢは、図２に示されている。これは、期間Ｔで周期的に発出される電磁／光パルス２０からなり（左側）、そのスペクトル２１は、周波数コムを含む（右側）。何れの周期性もＴＤＳ法の実行にとって重要である。

【0005】

今日、電磁（ＴＨｚ）パルスは、フェムト秒レーザパルスの光整流又はフェムト秒レーザパルスによって照明される光伝導性ＴＨｚアンテナによって取得され得る。

【0006】

ＴＨｚ範囲内にあるために、パルスは、数百フェムト秒～数ピコ秒にわたり継続すべきである。レーザの繰返し周波数ｆｒｅｐ（１／Ｔに等しい）は、ＧＨｚ～サブＨｚ周波数で様々であり得る。

【0007】

次に、レンズ又は放物面鏡を用いてパルスビームをコリメートし、それを調査対象のサンプルＳに向かって方向付ける。サンプルから発せられる透過（又は反射）パルスは、その後、光学システムによって収集されて、検出器Ｄに偏向される。検出器は、電磁波の電界を時間の関数としてフェムト秒～数百ピコ秒又はさらにナノ秒に及ぶスケールで測定する。これは、光伝導性又は電気光学サンプリングによって行われ得る。時間サンプリングは、典型的には、電動移動ステージ上に取り付けられたミラーで構成される遅延線ＤＬを用いてコヒーレント検出によって行われ、図３に示されるように、最大時間エクスカーションｔｍａｘを有する遅延Δｔが導入される。時間サンプリングは、その繰返し率がわずかに異なる２つの周波数コムのビート測定を通しても実行できる（ＨＡＳＳＯＰＳと呼ばれる技術）。２回の測定間の時間間隔は、サンプリング期間ｔｓと呼ばれ、周波数サンプリングｆｓ＝１／ｔｓである。下記では、時間トレースと呼ばれる、電界測定が行われる時間範囲をｔｍａｘと呼ぶ。

【0008】

典型的なサンプル期間ｔｓは、１０～５０ｆｓである。平均エネルギーではなく、ＴＨｚパルスの電界を直接測定できることにより、波形の位相及び振幅の両方を把握することができ、したがってサンプルの吸光係数及び屈折率に関する情報が提供される。本願では、測定サンプル時間トレースをＥｓ（ｔ）と呼び、サンプルがない状態で測定された参照時間トレースをＥｒｅｆ（ｔ）と呼ぶ。

【0009】

材料分析の場合、材料パラメータを取得する通常の方法は、記録されたパルス時間トレースＥｓ（ｔ）及びＥｒｅｆ（ｔ）（サンプルがある場合及びない場合）のフーリエ変換を行うことである。線形概算では、これらの２つのスペクトルの比は、複素透過係数と呼ばれ、以下のように記述することができる。

【数1】

【0010】

【数2】

及び

【数3】

は、それぞれ時間領域信号Ｅｓ（ｔ）及びＥｒｅｆ（ｔ）のフーリエ変換であり、

【数4】

は、複素屈折率であり、実部が遅延に対応し、虚部が材料中への吸光に対応し、ｄは、測定しなければならないサンプルの厚さであり、ωは、式ω＝２πｆによって周波数に関連付けられる角周波数である。項

【数5】

は、２つの空気／材料界面に関する通常入射でのフレネル係数の積であり、

【数6】

は、サンプル中のファブリペロー多重反射を考慮に入れた項である。

【0011】

式（１）は、いわゆる順問題を設定し、

【数7】

及び

【数8】

が分かっている場合、

【数9】

を得ることができる。実験により、

【数10】

及び

【数11】

が得られるため、実際の関心対象は、「逆問題」であり、すなわち、

【数12】

及び

【数13】

が分かれば、

【数14】

を特定することができる。

【0012】

一般に、サンプルの関心対象の物理パラメータを特定するための第一の方法１０は、以下の原理的なステップ：
－ サンプルのある場合及びない場合の実験を実行するステップ、すなわちＥｓ（ｔ）及びＥｒｅｆ（ｔ）を測定するステップ（ラクトースに関する例が図４に示されており、領域５は、減衰正弦曲線後の振動に対応する）、
－ 両方の時間トレース

【数15】

及び

【数16】

のｔｍａｘ（測定値がある場合）に等しい時間ウィンドウについて離散フーリエ変換を計算し（図４の信号に基づく図５を参照されたく、ピーク６は、減衰正弦曲線５のフーリエ変換の主成分に対応する）、及び測定伝達関数

【数17】

を特定するステップ、
－ サンプルの厚さを測定するステップ、
－ 屈折率の実部及び／又は虚部

【数18】

を、（勾配フリー又は準ニュートン等の最小化アルゴリズムを使用して）エラー関数を最小化することによって抽出するステップであって、エラー関数は、係数エラーδρ及びモデル化された透過係数

【数19】

と測定透過係数

【数20】

との間の位相エラーδψから定義され、抽出は、全ての単一の信号周波数について行われる（ラクトースに関する図６を参照されたく、ｄ＝９００μｍ±２０μｍである）ステップ、
－ 屈折率がドルーデ・ローレンツモデル等のモデルに追従することを前提として、屈折率をフィッティングして材料パラメータを取得するステップ
を含む。材料パラメータは、典型的には、スペクトル線の強度及び共振周波数並びに／又は線幅である（以下を参照されたい）。

【0013】

ステップｉｉ）に関して、離散フーリエ変換（ＤＦＴ）は、周期的離散関数のフーリエ変換を特定するための数学的ツールである。ＤＦＴの計算は、ＦＦＴアルゴリズム（高速フーリエ変換）によって最も頻繁に実行される。

【0014】

サンプルと相互作用する電磁パルスには、サンプルの屈折率の実部

【数21】

に関する遅延及び屈折率の虚部

【数22】

に関する吸光が生じる。図４の遅延Ｒは、

【数23】

に関する時間遅延であり、図４の振動５は、サンプルの吸光を波として見たものを示し、ピーク６は、

【数24】

に関する周波数領域の吸光ピークである。

【0015】

ステップ（ｖ）に関して、ドルーデ・ローレンツモデルは、サンプル中において、光が下降調波振動子内の電子と相互作用するという前提に基づく。その運動方程式は、時間領域の指数関数によって減衰する正弦曲線で表すことができ（図４の５を参照されたい）、これは、周波数ドメインのローレンツ曲線に対応する。吸光スペクトル線は、したがって、図７に示されるローレンツ曲線を辿る。このモデルは、特に、ガス分光法において圧力が数十マイクロバールより高くなり次第、使用される。

【0016】

ローレンツ曲線を辿る各スペクトル線は、３つのパラメータ：最大値Ｍ、線幅γ、共鳴周波数ｆ０によって特徴付けられ、それぞれ物理パラメータｐ１、ｐ２及びｐ３に対応する（この場合、ｉ＝１～３である）。

【0017】

この第一の方法は、確実な結果を示しているが、幾つかの欠点がある。

【0018】

第一に、サンプル厚さの精密な測定が必要である。第二に、結果が因果律（この問題ではクラマース・クローニッヒの関係式の形態をとる）に従わず、したがって誤った結果が導かれ得る。

【0019】

上述の方法の欠点の幾つかを克服するために、Ｐｅｒｅｔｔｉの刊行物と呼ぶ刊行物“ＴＨｚ－ＴＤＳ ｔｉｍｅ－ｔｒａｃｅ ａｎａｌｙｓｉｓ ｆｏｒ ｔｈｅ ｅｘｔｒａｃｔｉｏｎ ｏｆ ｍａｔｅｒｉａｌ ａｎｄ ｍｅｔａｍａｔｅｒｉａｌ ｐａｒａｍｅｔｅｒｓ”ｆｏｒｍ Ｐｅｒｅｔｔｉ ｅｔ ａｌｌ，ＩＥＥＥ ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｔｅｒａｈｅｒｔｚ ｓｃｉｅｎｃｅ ａｎｄ ｔｅｃｈｎｏｌｏｇｙ，ｖｏｌ９，ｎｏ．２，Ｍａｒｃｈ ２０１９”は、測定されたＴＨｚパルスの初期時間領域データ比較に基づく第二の方法を記載している。このような方法を実施するための理論上の方法は、
－ （ｉ）持続時間ｔｍａｘについての測定サンプル時間トレースＥｓ（ｔ）及び測定参照時間トレースＥｒｅｆ（ｔ）を取得するステップ、
－ （ｉｉ）サンプルが、参照パルスを、モデル化されたパルスにどのように変換するかを説明するパラメータ｛ｐ_ｉ｝の組に依存するモデルＥ_{ｍｏｄｅｌ}｛ｐ_ｉ｝（ｔ）を特定するステップ、
－ （ｉｉｉ）モデル化されたパルスと、測定された（サンプル）パルスとの間の差のＬ^２ノルム（差の二乗の合計の平方根）に等しい、最小化のための目的関数Ｏｂｊ｛ｐ_ｉ｝

【数25】

を特定するステップ、
－ （ｉｖ）目的関数を最小化するパラメータｐｉの値の組を特定するステップ
となるであろう。

【0020】

この方法は、サンプルの厚さｄが分かっていることを必要とせず、因果律に従っており、なぜなら、厚さがパラメータとして得られ、またエラー計算が時間領域内で行われるためである。

【0021】

残念ながら、サンプルが、参照パルスを、モデル化されたパルスにどのように変換するかを説明するモデリングは、時間領域において知られていない。既知のモデルの何れも、周波数ドメイン内でのサンプル挙動：

【数26】

を説明するものである。

【0022】

しかしながら、式（２）の公式化の実践的利点は、関数のノルムがそのフーリエ変換のノルムと同じであると明示するパーセバルの定理によって得られ、これは、目的関数Ｏｂｊ｛ｐ_ｉ｝が、以下の式：

【数27】

を用いて計算され得ることを意味し、式中、
ω_ｍａｘ＝π／ｄｔであり、ここで、ｄｔは、サンプル期間ｔｓに対応し、
ω_ｍｉｎ＝－π／ｄｔである。

【0023】

したがって、式（３）では、目的関数は、周波数領域でのサンプル挙動のモデリングと、測定サンプル時間トレースＥｓ（ｔ）の離散フーリエ変換との間の差の関数として表現される。

【0024】

これは、極めて好都合であり、時間及び周波数の両方の領域での目的関数の計算が可能となる。この方法は、したがって、周波数領域で直接実行することができ、繰返しごとのフーリエ変換の計算が不要となる。

【0025】

ＤＦＴを使用することにより、上述の２つの方法は、それらのデータが周期的であり、期間が記録の時間の長さに対応することを前提とする。その結果、これらは、２つのパルス間の時間距離Ｔが遅延線の最大エクスカーションｔｍａｘと正確に等しいことを前提とし、これは、一般的に誤りである。

【0026】

周波数領域で実行されるが、時間領域でのフィットに対応する、Ｐｅｒｅｔｔｉの刊行物に記載されている第二の方法２０は、図８に示されており、以下のステップを含む：
－ ステップＡ０：測定サンプル時間トレースＥｓ（ｔ）を取得するステップ、
－ ステップＢ０：測定参照時間トレースＥｒｅｆ（ｔ）を取得するステップ、
－ ステップＣ０：ｔｍａｘに等しい時間ウィンドウについてＥｓ（ｔ）の離散フーリエ変換

【数28】

を特定するステップ、
－ ステップＣ０’：ｔｍａｘに等しい時間ウィンドウについてＥｒｅｆ（ｔ）の離散フーリエ変換

【数29】

を特定するステップ。

【0027】

これらのステップは、第一の方法１０のステップ（ｉ）及び（ｉｉ）と同じである。
－ ステップＤ０：物理（所定の）パラメータｐｉに依存するサンプルの物理挙動モデルからサンプル周波数モデル

【数30】

を特定するステップ。このモデルは、周波数領域でのサンプルのインパルス応答のモデリングを行う。

【0028】

典型的には、モデル

【数31】

は、参照時間トレース

【数32】

のフーリエ変換に、式（１）の複素屈折率に関する伝達関数

【数33】

を乗じることを含む。

【数34】

－ ステップＥ０：物理パラメータｐｉの組に最適化アルゴリズムを適用するステップであって、
－ Ｅ１ 物理パラメータを初期化するステップ、
－ その後、測定サンプル時間トレース

【数35】

の離散フーリエ変換とサンプル周波数モデル

【数36】

との間の差から目的関数Ｏｂｊ｛ｐｉ｝を特定するステップＥ２を、前記目的関数を最小化するパラメータｐｉの値を得るまで繰り返し実行するステップ
を含むステップ。例えば、目的関数が所定の閾値より小さくなるまで繰り返す。

【0029】

例えば、目的関数は、以下の差の二乗（最小二乗アルゴリズム）から特定される。

【数37】

【0030】

方法は、計算を幾つかのステップで実行し、時間領域でのフィットを実現するが、周波数領域での簡単な計算を使用する。

【0031】

バルク中実材料に関して、Ｐｅｒｅｔｔｉの刊行物は、複数のドルーデ・ローレンツモデルを実装しており、これは、電子共振子（マトリクス振動、振動電荷等）の組としてサンプルの誘電率を定義し、以下の誘電率関数を導く。

【数38】

式中、ε_∞は、関心対象範囲と比較した高周波数での誘電率であり、ω_ｐは、プラズマ周波数であり、γ_ｐは、減衰率であり、ｋ_ｍａｘは、考慮される振動子の数であり、ω_ｏｋ、γ_ｋ及びΔε_ｋは、それぞれｋ番目の振動子の共鳴角周波数、減速率（線幅）及び強さ（最大、誘電率単位で表される）である。

【0032】

項

【数39】

は、自由キャリアの存在のためにルドルーデ項である。

【0033】

式（６）は、

【数40】

の式（１）に導入され、モデル

【数41】

は、式（５）によって定義される。

【0034】

したがって、自由キャリアのない１つの均一な層の場合（ドルーデ項が０に等しい）、伝播は、フレネル係数及び複数のローレンツ振動子モデルを通して３×ｋ＋１の物理パラメータｐｉの組を用いてモデル化されることになる。

【0035】

最適化アルゴリズムでは、したがって、各スペクトル線ｋについてパラメータｐ１ｋ_ｏｐｔ、ｐ２ｋ_ｏｐｔ、ｐ３ｋ_ｏｐｔの３つの値の組プラスε_∞の値が得られ、目的関数が最小化される。

【0036】

メタサーフェスに関する

【数42】

及び

【数43】

の他のモデルもＰｅｒｅｔｔｉの刊行物に記載されており、これは、時間領域結合モード理論（ＴＤＣＭＴ）と呼ばれる。

【0037】

上述の方法は、前述の調波振動子モデルを使用して波の時間トレースを直接モデル化し、この振動子（固有周波数、減衰持続時間及び電磁波との結合）のパラメータを、モデルがサンプル時間トレースとマッチするように調整することを提案している。これを行うために、パーセバルの定理により、最適化アルゴリズムが周波数領域について実行される。

【0038】

Ｐｅｒｅｔｔｉの刊行物に記載されている方法には、様々な利点がある：（ｉ）サンプル厚さの正確な測定が不要である。実際、炭化水素又は半導体ウェハ等の材料のサブマイクロメートルレベルの精度での正確なサンプル厚さを取得することは、困難である。したがって、このステップを回避することは、真の改善である。（ｉｉ）屈折率モデリングの問題は、全体として分析され、したがって通常の周波数につき２つの値と比較して、フィットのためのわずかなパラメータでよく、この方法は、ノイズの影響をはるかに受けにくい。これにより、屈折率に関して非常に高い精度に到達できる。（ｉｉｉ）残留フィットエラーは、振幅ユニット内にあるため、このエラーの解釈を明確に行うことができ、それは、実験のよりよい理解及び実装されるモデルの見落としの発見の可能性につながる。（ｉｖ）フィットは、時間領域で行われるため（最適化は、周波数領域で行われる）、強力な吸光があっても位相が失われず、追加のステップが不要である。（ｖ）最後に、これにより、ドルーデ・ローレンツモデルを用いた材料パラメータの精密、且つ確実、且つ一貫した取得だけでなく、メタ材料のそれらも可能となる。

【0039】

あらゆる古典的な方法と同様に、Ｐｅｒｅｔｔｉの刊行物に記載されている方法にも幾つかの欠点がある。主なものは、分解能であり、これは、最大時間遅延ｔｍａｘの値によって限定される。

【0040】

実際、パーセバルの定理により、最適化は、周波数領域内で行われる。時間ウィンドウｔｍａｘで実行される離散フーリエ変換

【数44】

及び

【数45】

の計算（ＦＦＴアルゴリズムによる）は、スペクトル分解能に直接的な影響を与え（以下を参照されたい）、これは、スペクトル線の幅を測定するか、又は２つの近接するスペクトル線を区別することが可能であることを意味する。

【0041】

スペクトル分解能δｆ、すなわちスペクトル線幅を測定するか、又は２つの隣接する線を区別する能力は、通常のいわゆるフーリエ不確実性限界によって限定され、これは、置換線の総遅延時間の逆周波数に対応する：
δｆ＝１／ｔｍａｘ。

【0042】

ｔｍａｘが約１ｎｓ（実現可能な最良値）であるとき、スペクトル分解能は、最良で１ＧＨｚの範囲である。

【0043】

１／ｔｍａｘのこのスペクトル限界は、ＤＦＴ計算から直接得られる。

【0044】

関数ｆの標準フーリエ変換は、

【数46】

である。

【0045】

間隔［ｎ，ｎ＋１［で一定の関数ｆでは、離散フーリエ変換は、

【数47】

である。

【0046】

周期的関数ｆで周期についてｙ１～ｙＮのＮ個の離散値をとる場合、ＤＦＴは、以下の式によって与えられる。

【数48】

【0047】

時間及び周波数は、２つの共役的な大きさであるため、ＤＦＴのスペクトル分解能は、ハイゼンベルクフーリエ基準によって限定される。ＤＦＴは、フーリエ変換限定分解能δｆを有し、簡潔に言えば、

【数49】

である。

【0048】

実際に、ＤＦＴ計算は、ｔｍａｘと正確に等しい周期性をとるため、ｔｍａｘ直後の時点に対応する情報は、０直後の時点に対応する情報に等しいと考えられる。これは、よく知られたスペクトルエイリアシングである。例えば、刊行物Ｘｕ，Ｊ．，Ｙｕａｎ，Ｔ．，Ｍｉｃｋａｎ，Ｓ．，＆Ｚｈａｎｇ，Ｘ．Ｃ．（２００３），“Ｌｉｍｉｔ ｏｆ ｓｐｅｃｔｒａｌ ｒｅｓｏｌｕｔｉｏｎ ｉｎ ｔｅｒａｈｅｒｔｚ ｔｉｍｅ－ｄｏｍａｉｎ ｓｐｅｃｔｒｏｓｃｏｐｙ”，Ｃｈｉｎｅｓｅ Ｐｈｙｓｉｃｓ Ｌｅｔｔｅｒｓ，２０（８），１２６６は、ＳＮＲ（信号対ノイズ比）が１より高いと仮定して、時間がｔｍａｘに限定された電磁パルスの検出が式（７）の分解能の周波数限定につながると説明している。

【先行技術文献】

【非特許文献】

【0049】

【文献】“ＴＨｚ－ＴＤＳ ｔｉｍｅ－ｔｒａｃｅ ａｎａｌｙｓｉｓ ｆｏｒ ｔｈｅ ｅｘｔｒａｃｔｉｏｎ ｏｆ ｍａｔｅｒｉａｌ ａｎｄ ｍｅｔａｍａｔｅｒｉａｌ ｐａｒａｍｅｔｅｒｓ”ｆｏｒｍ Ｐｅｒｅｔｔｉ ｅｔ ａｌｌ，ＩＥＥＥ ｔｒａｎｓａｃｔｉｏｎｓ ｏｎ Ｔｅｒａｈｅｒｔｚ ｓｃｉｅｎｃｅ ａｎｄ ｔｅｃｈｎｏｌｏｇｙ，ｖｏｌ９，ｎｏ．２，Ｍａｒｃｈ ２０１９”

【文献】Ｘｕ，Ｊ．，Ｙｕａｎ，Ｔ．，Ｍｉｃｋａｎ，Ｓ．，＆Ｚｈａｎｇ，Ｘ．Ｃ．（２００３），“Ｌｉｍｉｔ ｏｆ ｓｐｅｃｔｒａｌ ｒｅｓｏｌｕｔｉｏｎ ｉｎ ｔｅｒａｈｅｒｔｚ ｔｉｍｅ－ｄｏｍａｉｎ ｓｐｅｃｔｒｏｓｃｏｐｙ”，Ｃｈｉｎｅｓｅ Ｐｈｙｓｉｃｓ Ｌｅｔｔｅｒｓ，２０（８），１２６６

【発明の概要】

【発明が解決しようとする課題】

【0050】

したがって、サンプルの物理パラメータを、前述のスペクトル限界を克服し、したがって改善されたスペクトル分解能を有するＴＤＳ技術が実行されたそのサンプルから生じる電磁界の測定から特定する新規な方法が求められている。

【課題を解決するための手段】

【0051】

第一の態様によれば、サンプルの物理パラメータの組を特定する方法が提供され、これは、
－ Ａ 測定サンプル時間トレースＥｓ（ｔ）を取得するステップであって、測定サンプル時間トレースＥｓ（ｔ）は、時間領域分光法により、期間Ｔを有し、且つコム周波数を示す電磁パルスを周期的に発出する励起ビームによってサンプルを照明し、及びサンプルに由来する電磁場をコヒーレント検出によって時間の関数として検出することによって取得されており、サンプル時間トレースが測定される持続時間は、ｔｍａｘであり、ここで、ｔｍａｘ＜Ｔである、ステップ、
－ Ｂ 測定参照時間トレースＥｒｅｆ（ｔ）を取得するステップであって、測定参照時間トレースＥｒｅｆ（ｔ）は、ステップＡと同じ条件において、ただしサンプルの存在なしでの照明及び検出によって取得されている、ステップ、
－ Ｃ 期間Ｔにおいて広がり、且つ測定が行われなかった時点にゼロ値を付与することによって得られる、Ｅｒｅｆ０（ｔ）と呼ばれる拡張参照時間トレースを特定し、及びＴに等しい時間ウィンドウについて計算される、拡張参照時間トレースの離散フーリエ変換

【数50】

を特定するステップ、
－ Ｄ 拡張参照時間トレースのフーリエ変換

【数51】

及びサンプルの物理挙動モデルから、物理パラメータの組に依存する、サンプル周波数モデル

【数52】

と呼ばれる、周波数領域でのサンプルのインパルス応答のモデリングを特定するステップ、
－ Ｅ 物理パラメータの組に最適化アルゴリズムを適用するステップであって、
－ Ｅ１ 物理パラメータを初期化するサブステップ、
－ 以下のサブステップ：
－ Ｅ２ 推定サンプル時間トレースＥ_ｅｓｔ｛ｐ_ｉ｝（ｔ）と呼ばれる、サンプル周波数モデル

【数53】

の逆離散フーリエ変換を計算するサブステップ、
－ Ｅ３ 測定サンプル時間トレースＥｓ（ｔ）と推定時間トレースＥｅｓｔ（ｔ）との間の差からエラー関数を計算するサブステップ
を、前記エラー関数を最小化する物理パラメータの値の組を得るまで繰り返し実行するサブステップ
を含むステップ
を含む。

【0052】

第一の態様の発展形態によれば、励起ビームは、ＴＨｚ領域内にあり、１００ＧＨｚ～３０ＴＨｚに含まれる周波数を有する。

【0053】

第一の態様のさらなる発展形態によれば、最大時間遅延ｔｍａｘは、測定参照時間トレースのエネルギーの９５％超を含むように選択される。

【0054】

第一のもののさらなる発展形態によれば、サンプル周波数モデル

【数54】

は、フーリエ変換Ｅｒｅｆ０（ω）に、サンプル挙動を特徴付ける伝達関数Ｔ（ω）を乗じることを含む。

【0055】

好ましくは、伝達関数Ｔ（ω）は、複素屈折率ｎ（ω）に依存する。

【0056】

第一の態様のさらなる発展形態によれば、誘電率ε（ω）と呼ばれる複素屈折率の二乗は、各スペクトル線に関するドルーデ・ローレンツモデルに追従し、スペクトル線は、３つのパラメータ、振幅、減衰率と呼ばれる幅及び中心周波数の組によって特徴付けられる。

【0057】

第一の態様のさらなる発展形態によれば、エラー関数は、

【数55】

として定義される。

【0058】

第二の態様によれば、サンプルを特性評価するための特性評価装置が提供され、前記装置は、
－ 測定サンプル時間トレースＥｓ（ｔ）及び測定参照時間トレースＥｒｅｆ（ｔ）を記憶するメモリであって、測定サンプル時間トレースＥｓ（ｔ）は、時間領域分光法により、期間Ｔを有し、且つコム周波数を示す電磁パルスを周期的に発出する励起ビームによってサンプルを照明し、及びサンプルに由来する電磁場をコヒーレント検出によって時間の関数として検出することによって取得されており、サンプル時間トレースが測定される持続時間は、ｔｍａｘであり、ここで、ｔｍａｘ＜Ｔであり、
測定参照時間トレースＥｒｅｆ（ｔ）は、測定サンプル時間トレースＥｓ（ｔ）と同じ条件において、ただしサンプルの存在なしでの照明及び検出によって取得されている、メモリ、
－ 処理ユニットであって、
－ 期間Ｔにおいて広がり、且つ測定が行われなかった時点にゼロ値を付与することによって得られる、Ｅｒｅｆ０（ｔ）と呼ばれる拡張参照時間トレースを特定し、及びＴに等しい時間ウィンドウについて計算される、拡張参照時間トレースの離散フーリエ変換

【数56】

を特定すること、
－ 拡張参照時間トレースのフーリエ変換Ｅｒｅｆ０（ω）及びサンプルの物理挙動モデルから、物理パラメータの組に依存する、サンプル周波数モデル

【数57】

と呼ばれる、周波数領域でのサンプルのインパルス応答のモデリングを特定すること、
－ 物理パラメータに最適化アルゴリズムを適用することであって、
－ 物理パラメータを初期化するステップ、
－ 以下のサブステップ：
^＊推定サンプル時間トレースＥ_ｅｓｔ｛ｐ_ｉ｝（ｔ）と呼ばれる、サンプル周波数モデル

【数58】

の逆離散フーリエ変換を計算するサブステップ、
^＊測定サンプル時間トレースＥｓ（ｔ）と推定時間トレースＥｅｓｔ（ｔ）との間の差からエラー関数を計算するサブステップ
を、前記エラー関数を最小化する物理パラメータの値の組を得るまで繰り返し実行するステップ
を含む、適用すること
を行うように構成された処理ユニット
を含む。

【0059】

第三の態様によれば、分光光度計であって、
－ 本発明の第二の態様による特性評価装置、
－ 測定装置であって、
－ サンプルを励起ビームによって照明するように構成された光源、
－ 測定サンプル時間トレースＥｓ（ｔ）及び測定参照時間トレースＥｒｅｆ（ｔ）を検出するように構成された検出器
を含む測定装置
を含む分光光度計が提供される。

【0060】

以下では、下記の添付の図面を参照して本発明の実施形態並びにそのさらなる目的又は利点を詳細に説明する。

【図面の簡単な説明】

【0061】

【図1】ＴＨｚ時間領域分光法の原理を示す。

【図2】サンプルを照明する励起ビームの時間（左）及び周波数（右）の挙動を示す。

【図3】電磁パルスの時間サンプリングを示す。

【図4】時間の関数としてのサンプルの測定時間トレースＥｓ（ｔ）及び参照時間トレースＥｒｅｆ（ｔ）の例を示す。

【図5】時間トレースＥｓ（ｔ）及びＥｒｅｆ（ｔ）の両方の、ｔｍａｘに等しい時間ウィンドウについての離散フーリエ変換

【数59】

及び

【数60】

を示す。

【図6】エラー関数を最小化する屈折率の実部

【数61】

及び虚部

【数62】

を示す。

【図7】ローレンツ曲線を示す。

【図8】先行技術で説明されている方法を示す。

【図9】時間信号のＦＦＴによって得られた時間領域（左側）及びスペクトル領域（右側）の参照及びサンプル信号を示す。図９ａ（下）は、ｔｍａｘのスケールでの時間信号及びｔｍａｘに等しい時間ウィンドウで得られた周波数信号を示し、図９ｂ（上）は、各種のＴのスケールでの時間信号及びＴに等しい時間ウィンドウで得られた周波数信号を示す。

【図10】本発明による方法を示す。

【図11A】１ｋＰａの圧力でのサンプルの透過スペクトルを示す。

【図11B】約ｆ０＝０．５７３ＴＨｚの吸光スペクトルの拡大図を示す。

【図12】異なる方法：本発明による方法、従来技術、データベースによって得られた圧力Ｐの関数としての中心周波数の値の比較を示す。

【図13】異なる方法：本発明による方法、従来技術、データベースによって得られた圧力Ｐの関数としての減衰率γ（ＦＷＨＭとも呼ばれる）の値を示す。

【図14】本発明による、サンプルを特性評価するための特性評価装置５０を示す。

【図15】本発明による分光計を示す。

【図16】測定装置がコンピュータＩ／Ｏインタフェースに接続され、特性評価装置がコンピュータ内にある、本発明の実施形態を示す。

【図17】測定装置がＩ／Ｏインタフェースを介してコンピュータに接続され、特性評価装置が、コンピュータＣｏｍｐの通信サブシステムにインターネットを介して接続されたリモートサーバ内にある、他の実施形態を示す。

【発明を実施するための形態】

【0062】

本発明者らは、ＴＤＳ実験の詳細な分析を行い、従来技術の方式では常に時間ウィンドウｔｍａｘについてＤＦＴが実行されるという初期の所見を得た。その結果、前述のように、励起ビームＥＢの光パルスの周期性がｔｍａｘに等しいと仮定されることになる。

【0063】

この仮定は、実際の実験条件では誤りである。典型的には、２つのパルス間の期間Ｔは、ｔｍａｘより大きい。現在利用可能なレーザの場合、Ｔは、ｓ～ｎｓの範囲であり（ＧＨＺ～Ｈｚの繰返し率に対応する）、すなわち、ｔｍａｘは、典型的には、Ｔの２０％未満である。市販のレーザの繰返し率は、容易に変更することができず、非常に高い繰返し率、すなわちＴが１０ｎｓ未満であるレーザは、極めて高価であり、扱いにくい。

【0064】

したがって、特定の方法では、従来技術で行われているＤＦＴの計算は、経験と一致しない。その主な結果は、スペクトル限界が１／ｔｍａｘに等しいことである。

【0065】

分解能の改善は、図９に示されるように、測定が全期間Ｔについて行われた場合に得ることができる。左側は、時間領域の信号を示し、右側は、時間信号のＦＦＴによって得られるスペクトル領域の曲線を示す。

【0066】

信号Ｔｒは、シミュレーションによる参照時間トレースＥｒｅｆ（ｔ）の例であり、信号Ｔ１～Ｔ３は、Ｔ１からＴ３に減少する周波数幅Δｖを示す吸光線を有するサンプルで得られた、シミュレーションによるサンプル時間トレースＥｓ（ｔ）の例（すなわちＴ１からＴ３への時間領域のトレースの拡張）である。

【0067】

図の下側部分ａ）では、信号は、１０ｐｓの時間ｔｍａｘについてのみ入手でき、図の上側部分ｂ）では、信号は、期間Ｔ（１００ｐｓに等しい）より長い持続時間（３００ｐｓ）について入手できた。

【0068】

下側部分ａ）のスペクトルは、ｔｍａｘに等しい時間ウィンドウを用いてＦＦＴによって計算された。周波数領域の分解能Ｒｅｓは、限定され、Δｖの減少は、観察できないことが分かり、すなわち、スペクトル分解能は、１／ｔｍａｘ＝１００ＧＨｚに限定される。

【0069】

上側部分ｂ）のスペクトルは、Ｔに等しい時間ウィンドウを用いてＦＦＴによって計算された。周波数領域の分解能では、Δｖの減少を観察できることが分かり、すなわち、スペクトル分解能は、１／Ｔ＝１０ＧＨｚである。

【0070】

図の下側部分は、Ｔより短い持続時間ｔｍａｘについて行われた実際の実験の限界を示す。具体的には、この実験では、分解能は、時間の比Ｔ／ｔｍａｘに等しく、ここでは１０の係数で低下した。このスペクトル限界は、気体の狭い線の特性評価にとって非常に大きい欠点である。

【0071】

本発明による、サンプルの物理パラメータの組を特定する方法３０が図１０に示されている。方法３０は、測定サンプル時間トレースＥｓ（ｔ）を取得するステップＡと、測定参照時間トレースＥｒｅｆ（ｔ）を取得するステップＢとを含む。好ましくは、両方の測定時間トレース、サンプルの時間トレース及び参照時間トレース（ステップＡと同じ条件において、ただしサンプルの存在なしでの照明及び検出によって得られる）は、前述のように古典的なＴＤＳ装置によって得られる。測定サンプル時間トレースＥｓ（ｔ）は、周期Ｔの、周波数コムを示す電磁パルス（波長が１０μｍ未満である場合、光パルスと呼ばれ得る）を周期的に発出する励起ビームＥＢによってサンプルＳを照明して特徴付け、且つサンプルに由来する電磁場（透過又は反射）をコヒーレント検出によって時間の関数として検出することによって得られている。典型的には、コヒーレント検出は、遅延線を含むか、又はＨａｓｓｏｐｓ実験（デュアル周波数コム）に基づく。

【0072】

好ましくは、励起ビームＥＢは、ＴＨｚ領域にあり、０．１ＴＨｚ～３０ＴＨｚに含まれる周波数を有するが、特許請求される方法３０は、照明及び検出装置が利用可能となる場合、より高い周波数（ＩＲ、可視）も外挿法によって推定できる。

【0073】

両方の時間トレースが測定される持続時間は、ｔｍａｘと呼ばれ、ここで、ｔｍａｘは、Ｔより短い。

【0074】

特許請求される方法３０のステップＡ及びＢは、前述のＰｅｒｅｔｔｉの刊行物の方法２０のステップＡ０及びＢ０と同じである。

【0075】

ステップＣでは、Ｅｒｅｆ０（ｔ）と呼ばれる拡張された（「パディングされた」ともいわれる）参照時間トレースが特定される。Ｅｒｅｆ０（ｔ）は、期間Ｔにおいて広がり、且つ測定が行われなかった時点、すなわちｔ＝ｔｍａｘ後でｔ＝Ｔまでの時点にゼロ値を付与することによって得られる。これらの時点には、何れの任意の値も付与できるが、ゼロは、間隔］ｔｍａｘ，Ｔ］でほとんど又は全くパルスエネルギーが見られないことを表現するための最善の値である。この理由から、好ましくは、最大時間遅延ｔｍａｘは、測定参照時間トレースのエネルギーの９５％超（好ましくは９９％）を含むのに十分に大きくなるように選択される（方法の精度は、このパーセンテージの値と共に向上する）。データを計算する前に未知の値を０とするこの方法は、「０パディング」として知られる。

【0076】

また、ステップＣでは、拡張参照時間トレースＥｒｅｆ０（ｔ）の離散フーリエ変換

【数63】

が特定され、計算は、Ｔに等しい時間ウィンドウについて実行され、これは、Ｅｒｅｆ０（ｔ）が期間Ｔ全体にわたることから可能である。

【数64】

【0077】

したがって、

【数65】

を計算するＤＦＴは、時間ウィンドウｔｍａｘについて実行される、方法２０の、

【数66】

を計算するステップＣ０’のＤＦＴと異なる。

【0078】

次に、ステップＤでは、サンプル周波数モデル

【数67】

と呼ばれる、周波数領域でのサンプルのインパルス応答のモデリングが特定される。サンプル周波数モデルは、物理パラメータｐｉ（ｉは、パラメータのインデックスである）の組に依存し、拡張参照時間トレース及びサンプルの物理挙動モデルのフーリエ変換

【数68】

から得られる。ステップＤのモデリングは、方法２０のステップＤ０のモデリングと同じであるが、特許請求される方法３０のモデリングは、

【数69】

と異なる

【数70】

から定義される点が異なる。

【0079】

好ましくは、サンプル周波数モデル

【数71】

は、フーリエ変換

【数72】

に、サンプル挙動を特徴付ける伝達関数Ｔ（ω）を乗じることを含む。

【数73】

【0080】

好ましくは、伝達関数Ｔ（ω）は、前述のように、複素屈折率ｎ（ω）に依存する（式（１）を参照されたい）。

【0081】

モデルの選択は、サンプルに依存するが、分光法による線の形状が物理の法則によって決められるため、制約される。

【0082】

１つの実施形態によれば、誘電率ε（ω）と呼ばれる複素屈折率の二乗は、各スペクトル線に関するドルーデ・ローレンツモデルに追従し（ただし、それに限定されない）、すなわち、各スペクトル線は、ローレンツ分布の形状に追従する。この場合、スペクトル線ｋは、式（６）で説明されているように、３つのパラメータ：振幅ｐ１＝Δε_ｋ（誘電単位における）、減衰率と呼ばれる幅ｐ２＝γ_ｋ及び中心（共鳴）周波数ｐ３＝ω_０ｋの組によって特徴付けられる。このモデルは、気体のスペクトル線のモデリングに適している。

【0083】

ステップＥでは、最適化アルゴリズムが物理パラメータｐｉの組に適用される。まず、サブステップＥ１では、物理パラメータｐｉが初期化され、その後、以下のサブステップＥ２及びＥ３が繰り返し実行される。

【0084】

Ｅ２では、サンプル周波数モデルの逆離散フーリエ変換

【数74】

が計算される。この逆離散フーリエ変換は、推定サンプル時間トレースＥ_ｅｓｔ｛ｐ_ｉ｝（ｔ）と呼ばれ、以下のように特定される。

【数75】

ｔｓは、サンプル時間である。

【0085】

ここで、実験の実際の期間、すなわちレーザ源の繰返し率Ｔが計算に使用される。そのため、モデル及び実験は、同じ周期性に追従し、あらゆるＤＦＴ折り返し（エイリアシング）アーチファクトを防止するが、実際のものを再現する。その結果、フーリエ－ハイゼンベルク限界より狭い線では、モデリングによって導入される追加情報により、時間信号が時間ウィンドウの右側の縁を越え、左側の縁でその最初に戻ることになる。依然として、この信号は、フィットされ、そのため、フーリエ－ハイゼンベルク基準より高い分解能が得られる。

【0086】

次に、Ｅ３では、エラー関数εｅｒ｛ｐｉ｝が測定サンプル時間トレースＥｓ（ｔ）と推定時間トレースＥｅｓｔ（ｔ）との間の差から計算され、すなわち、εｅｒ｛ｐｉ｝は、関数ｆに従って［Ｅｓ（ｔ）－Ｅｅｓｔ（ｔ）］に依存し、ｆは、形態学的距離を定義する関数である。

【数76】

【0087】

繰返しは、エラー関数を最小化するパラメータｐｉの値の組を得るまで続く。例えば、繰返しは、エラー関数が所定の閾値より小さくなるか、又は展開が止まったときに停止する。

【0088】

終了した最適化の結果は、パラメータｐｉの値の組である。

【0089】

ローレンツモデルに追従する気体スペクトル線の場合、線ごとのパラメータの３つの値の組（Δε_ｋ、γ_ｋ、ω_０ｋ）が提供される。ｐ３＝ω_０ｋにより、気体の識別が可能となり、ｐ１＝Δε_ｋは、気体の圧力に関係し、ｐ２＝γ_ｋは、気体の温度に関係する。

【0090】

最適化が０～ｔｍａｘのｔの値について行われる点に留意することが重要である。記録されたデータＥｓ（ｔ）は、ｔｍａｘについてのみ入手可能である（時間フレームＴ全体についての測定は、前述のように不可能である）。Ｅｓ（ｔ）上に存在しない情報を追加することができないため、モデルの時間フレームの一部分のみが、記録されたデータＥｓ（ｔ）と比較され得る。典型的には、Ｅｅｓｔ（ｔ）の利用可能な算出データの約１０％（ｔｍａｘ／Ｔの比）がエラー計算に使用される。

【0091】

その結果、仮定が満たされないため、Ｐｅｒｅｔｔｉの刊行物の方法２０で以前に行われていたように、パーセバルの定理を用いて曲線を比較し、周波数領域のエラーを計算することができない。

【0092】

パーセバルの定理を使用できないため、ステップ２における各繰返しでは、時間領域のモデル化されたデータの取得及びその後のエラー計算には、高速フーリエ変換を行わなければならない（ステップＥ３）。

【0093】

エラー関数の選択は、最適化の種類に依存する。エラー関数の例は、二乗平均の平方根の差である。

【数77】

【0094】

例として、拡張ラグランジュ粒子群オプティマイザが条件付き最適化を実行するために使用される。

【0095】

サンプル挙動、すなわち、

【数78】

が追従するか、又は該当する場合にはＴ（ω）が追従する式は、既知である。この「アプリオリ」知識は、非常に重要であり、なぜなら、それにより、この条件付き再構成アルゴリズムを何れの情報も失わずに実行できるからであり、それにより、分解能は、Ｐｅｒｅｔｔｉの刊行物の方法２０よりはるかに高くなる（超分解能）。

【0096】

例えば、気体のスペクトル線の理論物理学によって与えられる知識から、記録される線は、時間領域の指数関数的減衰正弦曲線に対応する周波数領域内のローレンツ分布形状を辿ると言うことができる。超分解能は、各減衰正弦曲線が３つのパラメータ：振幅、中心周波数及び減衰率のみによって説明されることに依拠する。パラメータを取得するためにわずかな時間領域点があればよく、全時間トレースの点の総数の一部より振動子の数が少ないという限定がある（一般化フーリエ－ハイゼンベルク不確実性）。換言すれば、ＴＤＳシステムは、情報として全時間トレースを記録し、それから３ｋ＋１のパラメータを取得することが多くの場合に望まれる。時間領域内で計算を行うことは、周波数領域内で各δｆ（δｆは、フーリエ変換の周波数ステップである）の情報の検索を行う代わりに、わずかなモード（ｋ個の振動子）が検索され、したがってスパース情報とも呼ばれるはるかに少ない情報が検索されることを意味する。この制約により、分解能は、超分解能微鏡検査法と同様に、信号対ノイズ比及びフィットを最適化する計算方法によってのみ限定されることになる。

【0097】

特許請求されるサンプルパラメータ特定方法３０は、非常に広いスペクトルについて高分解能分光法を実行することにより、ＴＤＳスペクトル内に存在する情報量を十分に利用することを可能にする。それにより、例えば気体の数百の線（同じ気体又は混合物）を検査することにより、したがって異なる成分の相対濃度を測定することによりイベントをモニタすることができる。

【0098】

多くの計算時間（１回の繰返しにつき１回のＤＦＴ）を費やす、特許請求される方法の速度を高めるために、特許請求される方法の１つの実施形態において、制約アルゴリズムの代わりに高調波反転方式がとられる。高調波反転は、全ての要件が満たされた場合に使用され得、特に励起パルスがディラックのパルスにごく近くなるようにして、システムを時間トレースの大部分で十分に弛緩させることが確実となるようにすることが重要である。

【0099】

方法のパフォーマンスを試験するために、アンモニア（ＮＨ_３）の時間トレースを、小型の可搬式オールファイバ型分光計である市販のＴＨｚ－ＴＤＳ ＴＥＲＡＳＭＡＲＴ Ｍｅｎｌｏ（登録商標）によって記録した。これは、繰返し率１００ＭＨＺのフェムト秒レーザ（９０ｆｓパルス）を使用する。ＴＨｚパルスの持続時間は、２００ＧＨｚ～５ＴＨｚのスペクトルで約４００ｆｓである。光路長８ｃｍのブリュースタアングルシリコンウィンドウガスセルを光路内にセットして測定を行った。１３分間の実験を３ｍｂ～１００ｍｂの異なる圧力で繰り返した。１ｋＰａ（１０ｍｂ）でのサンプルの透過スペクトルで、図１１Ａに示されるスペクトルが得られた。主な吸光ピークは、ｆ０＝０．５７３ＴＨｚでのピークである。

【0100】

図１１Ｂは、ｆ０＝０．５７３ＴＨｚのピーク及びその前後のＦＦＴの拡大図を示す。周波数サンプリングｆｓが１．２ＧＨＺに等しい状態では、ピークは、１点でのみ見えることが分かる。

【0101】

図１２は、特許請求される方法３０によって得られる、圧力Ｐに関する中心周波数の値ｆ０を示し（曲線１２１、三角）、図１３は、特許請求される方法３０によって得られる、圧力Ｐに関する減衰率γ（ＦＷＨＭとも呼ばれる）の値を示す（曲線１３１、三角）。

【0102】

データベースのＨｉｔｒａｎ ２０１６から抽出されたＮＨ_３のデータも比較のために図１２及び１３にプロットされている（それぞれ曲線１２２及び１３２、四角）。方法３０によって特定された物理パラメータの値ｆ０及びγは、γ及びｆ０の何れについても、Ｈｉｔｒａｎからのデータに非常によく適合し、したがって特許請求される方法のきわめて高い精度を実証している。

【0103】

Ｐｅｒｅｔｔｉの方法２０で得られたｆ０及び減衰率の値も比較のためにプロットされている（それぞれ曲線１２３及び１３３、丸）。図１３では、方法２０によってＦＷＨＭの高い値で良好な結果が得られることが分かり、これは、高い圧力に対応するが、１０ｍｂ未満の圧力でのＦＷＨＭの値は、方法２０の不十分な周波数分解能のため、格段に精度が低くなる。

【0104】

ＮＨ_３の約５３０ＧＨｚに焦点を当てることに加えて、スペクトルの他の線の分析を行い、５３０ＧＨｚ～３ＴＨｚ超の全てのスペクトルに沿った超分解能を得た。より具体的には、約１２５０ＧＨｚの線に焦点を当てた。この線は、実際に、３５０ＭＨｚ（遅延線の「フーリエ変換」分解能限界の３分の１）だけ分離されるダブレットである。特許請求される超分解能の方法により、ダブレットの各線の周波数は、数ｍＢａｒまでの圧力について、３０ＭＨｚより高い精度で得られている。したがって、分解能限界より低いダブレット分割が実現されており、これは、超分解能であることの確固たる証拠である。

【0105】

特許請求される方法の高い分解能から、ＴＨｚパルスへのＴＤＳ超分解能方法３０の適用により、ＴＨｚ－ＴＤＳの用途が特に呼吸アナライザのような大気又は健康の目的のため及び産業環境の管理のための気体分光法に広がる。

【0106】

他の態様によれば、本発明は、図１４に示される、サンプルＳを特性評価するための特性評価装置５０に関する。装置５０は、既に定義した測定サンプル時間トレースＥｓ（ｔ）及び測定参照時間トレースＥｒｅｆ（ｔ）を記憶するメモリＭＥＭと、処理ユニットＰＵであって、
－ 期間Ｔにおいて広がり、且つ測定が行われなかった時点にゼロ値を付与することによって得られる、Ｅｒｅｆ０（ｔ）と呼ばれる拡張参照時間トレースを特定し、及びＴに等しい時間ウィンドウについて計算される、拡張参照時間トレースの離散フーリエ変換

【数79】

を特定すること、
－ 物理パラメータｐｉの組（ω）に依存する、サンプル周波数モデル

【数80】

と呼ばれる、周波数領域でのサンプルのインパルス応答のモデリングを拡張参照時間トレースのフーリエ変換

【数81】

及びサンプルの物理挙動モデルから特定すること、
－ 物理パラメータ（ｐｉ）の組に対して最適化アルゴリズムを適用することであって、
－ 物理パラメータ（ｐｉ）を初期化するサブステップ、
－ 以下のサブステップ：
－ 推定サンプル時間トレースＥ_ｅｓｔ｛ｐ_ｉ｝（ｔ）と呼ばれる、サンプル周波数モデル

【数82】

の逆離散フーリエ変換を計算するサブステップ、
－ 測定サンプル時間トレースＥｓ（ｔ）と推定時間トレースＥｅｓｔ（ｔ）との間の差からエラー関数（ε_ｅｒ｛ｐｉ｝）を計算するサブステップ
を、前記エラー関数を最小化する物理パラメータの値（ｐｉ_ｏｐｔ）の組を得るまで繰り返し実行するサブステップ
を含む、適用すること
を行うように構成された処理ユニットＰＵとを含む。

【0107】

他の態様によれば、本発明は、図１５に示される分光計Ｓｐｅｃｔｒｏに関し、これは、ｉ）サンプルＳを励起ビームＥＢによって照明するように構成された光源ＬＳと、測定サンプル時間トレースＥｓ（ｔ）及び測定参照時間トレースＥｒｅｆ（ｔ）を検出するように構成された検出器Ｄとを含む測定装置ＭｅＤと、ｉｉ）特性評価装置５０とを含む。

【0108】

上述の実施形態は、単に非限定的な例であると理解されたい。特に、測定装置ＭｅＤ及び特性評価装置５０は、異なる要素内に配置され、何れの組合せでも一緒に使用され得る。

【0109】

ある実施形態において、図１６に示されるように、ＭｅＤは、Ｉ／Ｏインタフェース７０３を介してコンピュータＣｏｍｐに接続され得、特性評価装置５０ ＤＢＰは、コンピュータ内に配置される。他の実施形態では、図１７に示されるように、ＭｅＤは、Ｉ／Ｏインタフェース７０３を介してコンピュータＣｏｍｐに接続され得、特性評価装置５０は、コンピュータＣｏｍｐの通信サブシステム７２０にインターネット７５を介して接続されたリモートサーバ７６に配置され得る。

【0110】

他の態様において、本発明は、特許請求される方法のステップを実施するように適合されたコンピュータプログラムに関する。他の態様では、本発明は、そのコンピュータプログラムを組み込んだコンピュータ可読媒体に関する。

【図1】

【図2】

【図3】

【図4】

【図5】

【図6】

【図7】

【図8】

【図9】

【図10】

【図11A】

【図11B】

【図12】

【図13】

【図14】

【図15】

【図16】

【図17】