超伝導遷移端検出器 TES の I-V 解析：基礎方程式から実験方法まで

１ はじめに

　超伝導転移端センサ（Transition Edge Sensor; TES） の電流電圧（I-V）特性を測定することは、デバイスの基本性能を評価する上で欠かせないプロセスです。I-Vデータ解析を通じて、TESの物理的な特性を数値として明らかにすることができます。
　本記事では、TESのI-Vデータから諸パラメータを決定するための実践的な解析フローについて、初学者向けに詳しく解説します。

２ 基本用語

• 超伝導転移端センサ（Transition Edge Sensor; TES）
  　超伝導から常伝導の相転移点近傍における、急峻な抵抗変化を利用した高感度な温度センサです。TESは入射したエネルギーを熱として受け取ることで、その熱による温度上昇を抵抗値の変化に変換し、そのときに回路を流れる電流の変化を検出します。超伝導転移近傍の物理を利用するため、100mK程度の極低温環境で動作させる必要があります。
• 超電導量子干渉計（Sperconducting Quantum Interference Device; SQUID）
  　TESの性能を最大限に引き出すためには、読み出し回路にも高い性能が求められます。TESから出力される信号を損なうことなく電流変化を読み出すには、極低温環境下で動作することはもちろん、回路電圧への影響を最小限に抑える低インピーダンス特性と、微小な信号を埋もれさせない低ノイズ特性を備えた電流計が不可欠です。
  　これらの条件を満たすのが、SQUIDを用いた読み出し系です。SQUIDは微小な磁束変化を読み取る磁気センサですが、隣に磁束を入力するためのコイルを置くことによって、SQUIDを感度の高い電流計として扱うことができます。
• 電熱フィードバック（Electro-Thermal Feedback; ETF）
  　TESは温度計として高い感度を持っていますが、温度域が数mKと非常に狭いため、動作温度を超伝導転移端中に保つ必要があります。これは、電熱フィードバックと呼ばれる、TESを 定電圧バイアス で動作させ、強いフィードバックをかける技術で実現されます。

\rm 温度上昇 \longrightarrow 抵抗増加 \longrightarrow 電流減少（Ohmの法則） \longrightarrow 発熱減少

• Python
  　プログラミング言語の一つです。Pythonは文法がシンプルで読み書きしやすく、高い汎用性を有しているため広く利用されています。特に、実用的なライブラリが豊富に揃っている点が特徴です。今回のI-Vデータ解析もPythonを用いて行います。なお、実行するPython環境にはGoogle Colabを使用します。

３ 基礎知識

　いきなり解析内容の説明をする前に、I-Vデータ解析から得られるTESの諸パラメータの物理的意味を明らかにするとともに、評価に用いる特性グラフの評価、および測定の具体的な手順について順を追って説明します。

３.１ 温度計の感度

　２章で説明した通り、TESは超伝導転移端における急峻な抵抗変化を利用しています。温度計の性能を示す指標として、無次元量の感度 $\alpha$ を次のように定義します：

\alpha \equiv \frac{{\rm d}\log{R}}{{\rm d}\log{T}} = \frac{T}{R}\frac{{\rm d}R}{{\rm d}T}\ . \tag{1}

　ここで、 $T$ は温度計の温度、 $R$ は温度計の抵抗値です。この $\alpha$ は、温度変化に対して抵抗がどれだけ敏感に反応するかを表します。これが大きいほど、微小な温度変化、つまり外部からの微小なエネルギー入射を大きな抵抗値変化として変換できます。すなわち、ごくわずかなエネルギー量の差異を明確な電気信号の差として読みだせ、優れたエネルギー分解能を有していると言えます。
　TESの温度抵抗特性を以下に示します。このグラフが示すように、超伝導転移端における急激な抵抗値上昇で $\alpha$ が飛躍的に高められ、極めて高いエネルギー分解能が実現されます。

３.２ フィードバック

　まずは、TESの読み出し回路について説明します。いきなりですが、以下にその模式図を示します。少し複雑な回路に見えるかもしれませんが、大きく二つに分けて見ると理解しやすくなります。

　一つは、TESの動作温度を超電動転移端に安定して保ちつつ、外部からのエネルギー入力に伴う抵抗変化を電流信号として読み出すための、 シャント抵抗 を用いた擬似定電圧バイアス回路です。シャント抵抗をTESと並列に接続することで擬似的な定電圧バイアス環境が形成され、２章で説明した電熱フィードバックが実現されます。
　外部からの入力によってTESの温度が上昇すると、抵抗値の増加に伴って流れる電流が減少します。この電流変化は入力コイルを介して後段のSQUIDへ伝達され、検出信号として利用されます。また、抵抗値の増加によってJoule発熱が抑制されるため、電熱フィードバックが働き、TESは超伝導転移端に自動的に保持されます。

　二つは、前段の電流変化を高感度に捉え、電気信号として読み出すための、SQUIDの 磁束固定ループ（Flux Locked Loop; FLL） 回路です。入力コイルによって生じる磁束をSQUID素子で検出し、増幅された電圧信号として出力することで、２章で説明したSQUIDを用いた読み出し系が実現されます。
　SQUIDは外部磁束に対して周期的な応答をするため、動作点がわずかにずれるだけでも増幅率が大きく変化し、応答が非線形になります。そこでFLLでは、SQUIDに加わる磁束の変化を打ち消す向きのフィードバック磁束を印加し、SQUIDを貫く磁束が一定に保たれるようにします。

　フィードバック回路について少し説明しておきましょう。以下に示す図は、一般的なフィードバック回路です。直流増幅率 $A$ 、フィードバック量 $b$ とする増幅器を考えています。このとき、入力 $V_{\rm in}$ に対しての出力 $V_{\rm out}$ は

\begin{equation}
 \begin{aligned}
  (V_{\rm in} - bV_{\rm out})\times A=V_{\rm out}\ , \\
  \therefore \quad V_{\rm out} = \frac{1}{b}\frac{\mathcal L}{1+\mathcal L}V_{\rm in}\ ,\qquad \mathcal L = bA 
 \end{aligned}\tag{1}
\end{equation}

と計算できます。 $\mathcal L$ はループゲインと呼ばれ、特に $\mathcal L\gg 1$ のとき $V_{\rm out}\simeq V_{\rm in}/b$ となり、出力は $1/b$ 倍となることがわかります。

　これをSQUIDの磁束固定ループに適用します。SQUIDの出力 $V_{\rm out}$ は、フィードバック抵抗を介してフィードバックコイルに戻されるので、フィードバック量 $b$ はフィードバックコイルの作る磁束 $\Phi_{\rm FB}$ から $bV_{\rm out}=\Phi_{\rm FB}$ と書けます。さらに、SQUIDとフィードバックコイルの相互インダクタンス $M_{\rm FB}$ とフィードバック抵抗の抵抗値 $R_{\rm FB}$ を用いれば

b = \frac{\Phi_{\rm FB}}{V_{\rm out}} = \frac{M_{\rm FB}\cdot(V_{\rm out}/R_{\rm FB})}{V_{\rm out}} = \frac{M_{\rm FB}}{R_{\rm FB}} \tag{2}

と書き直せます。つまり、 $\mathcal L\gg1$ でのFLL回路のゲインは

\frac{1}{b} = \frac{R_{\rm FB}}{M_{\rm FB}} \tag{3}

です。
　SQUIDへの入力電流 $I_{\rm in}$ と出力電圧 $V_{\rm out}$ の関係については、 $|\Phi_{\rm in}|=|\Phi_{\rm FB}|$ より、式$(2)$から

\begin{align}
 \frac{M_{\rm in}I_{\rm in}}{V_{\rm out}} &= \frac{M_{\rm FB}}{R_{\rm FB}} \notag\\
 V_{\rm out} &= \frac{M_{\rm in}}{M_{\rm FB}}R_{\rm FB}I_{\rm in} \tag{4}
\end{align}

が得られます。このときの

\frac{M_{\rm in}}{M_{\rm FB}}R_{\rm FB} =: \Xi \tag{5}

を電流電圧変換係数と呼びます。

相互インダクタンス
　あるコイルに流した電流が、別のコイルにどれだけ磁束を結合させるかを表す比例定数です。コイル１に電流 $I_1$ を流すと磁場が発生し、その磁場の一部が近くに置かれたコイル２を貫きます。このとき、コイル２を貫く磁束 $\Phi_2$ は $\Phi_2\propto I_1$ になるので、比例定数をコイル１とコイル２の相互インダクタンス $M$ と定義して

\Phi_2 = MI_1 \tag{6}

と書きます。

３.３ I-V特性

　前節で述べたTESの読み出し回路の構成に基づき、I-V測定の理論的な取り扱いについて説明します。以下の記述は、３.２章で示した図中の文字をそのまま使います。また、TESに流れる電流を $I_{\rm TES}$ 、TESの両端にかかる電圧を $V_{\rm TES}$ で表記することとします。

３.３.１ 理論式

　I-V測定は、熱浴温度 $T_{\rm bath}$ を一定に保った状態で、Carolimeter Biasのバイアス電流 $I_{\rm bias}$ を走査させたときのSQUID出力電圧 $V_{\rm out}$ を取得することで行います。TESをシャント抵抗によ擬似定電圧バイアス下で動作させる場合、シャント抵抗とTESにかかる電圧は同じなので

\begin{align}
 R_{\rm sh}(I_{\rm bias} - I_{\rm TES}) &= R_{\rm TES}I_{\rm TES} \notag\\
 R_{\rm TES} &= \left(\frac{I_{\rm bias}}{I_{\rm TES}}-1\right)R_{\rm sh} \tag{7}
\end{align}

が成り立ちます。さらに、式$(4)$より出力電圧 $V_{\rm out}$ とTESに流れる電流 $I_{\rm TES}$ には

V_{\rm out}\simeq \Xi I_{\rm TES}\ ,\qquad \Xi =\frac{M_{\rm in}}{M_{\rm FB}}R_{\rm FB} \tag{8}

の関係があるので、以下の関係式が得られます：

R_{\rm TES}=\left(\Xi\frac{I_{\rm bias}}{V_{\rm out}}-1\right)R_{\rm sh}\ . \tag{9}

３.３.２ グラフの概形

　３.１章で見たように、TESの抵抗値は転移温度 $T_{\rm c}$ の前後で大きく特徴が異なります。ただ、TESに流れる電流とTESの両端にかかる電圧は、変わらずOhmの法則で書くことができます：

I_{\rm TES} = \frac{1}{R_{\rm TES}}V_{\rm TES}\ . \tag{10}

ゆえに、超伝導領域（$T<T_{\rm c}$）では、抵抗がゼロになるためI-V特性は傾きが1より十分大きい一次関数になることがわかります。また常伝導領域（$T>T_{\rm c}$）では、抵抗がTES固有の一定値になるため、I-V特性はある傾きを持った一次関数になります。
　転移温度付近（$T\approx T_{\rm TES}$）では、抵抗値が温度に依存して変化するため、I-V特性が一次関数にはなりません。しかし、転移温度付近では温度変化がほとんどないため、式$(6)$に注目すれば $P_{\rm TES}=I_{\rm TES}V_{\rm TES}$ が一定値になることがわかります。したがって、 $I_{\rm TES}$ は $V_{\rm TES}$ に反比例します。
　以上より、TESのI-V特性グラフは以下のような概形となります。

３.４ 熱伝導度

　I-V測定を行うことで得られる、TESの熱伝導度 $G$ とその温度依存性のべき定数 $n$ について説明します。ここで扱う熱伝導度とは、TESから熱浴へと流れる熱の伝わりやすさを指しており、文脈によっては「熱コンダクタンス」とも呼ばれます。
　極低温において、熱伝導率 $\kappa$ は温度 $T$ に対してべき乗則 $\kappa\propto T^{n-1}$ に従います（後の式変形を考慮して $n-1$ 乗としています）。例えば、電子が熱伝導を担う場合は $n=2$ 、phononが熱伝導を担う場合は $n=4$ です。熱伝導度は、熱が伝わる物質の長さを $\ell$ 、断面積を $A$ として

G = \frac{\kappa A}{\ell} \tag{11}

と書けますから、比例定数 $G_0$ で

G = G_0T^{n-1} \tag{12}

と表せます。ゆえに発熱量 $P$ と温度 $T$ の関係は非線形となります。そのため、熱伝導度を熱抵抗の逆数として単純に $G=P/T$ とするのではなく

G \equiv \frac{{\rm d}P}{{\rm d}T} \tag{13}

で定義します。
　一般に、熱浴の温度 $T_{\rm bath}$ に対して $T\gg T_{\rm bath}$ なので、TESの熱は熱浴へと流れ、このときに流れる熱量 $P_{\rm TES-nath}$ は

\begin{align}
 P_{\rm TES-bath} &= \int_{T_{\rm bath}}^{T} G\ {\rm d}T \notag\\
 &= \int_{T_{\rm bath}}^{T} G_0T^{n-1}\ {\rm d}T \notag\\
 &= \frac{G_0}{n}(T^n - T_{\rm bath}^n) \tag{14}
\end{align}

と計算できます。TESの温度が $T_{\rm TES}$ で一定となり、一定の熱量が熱浴へ流れる平衡状態になると、TESのJoule発熱 $P_{\rm TES}=V_{\rm TES}I_{\rm TES}$ とつり合うので

P_{\rm TES} = V_{\rm TES}I_{\rm TES} = \frac{G_0}{n}(T_{\rm TES}^n - T_{\rm bath}^n) \tag{15}

が得られます。ただし、 $V_{\rm bias}$ はバイアス電圧、 $R_{\rm TES}$ は動作点でのTESの抵抗値です。また、 $T_{\rm TES}$ について解くことで

T_{\rm TES} = \left(T_{\rm bath}^n + n\frac{P_{\rm TES}}{G}\right)^{\frac{1}{n}} \tag{16}

も得られます。

３.５ 測定方法

　最後に、実際に$I-V$測定を行う際の手順について簡単に説明します。

1. SQUIDの$\phi$-Vが最大となる値にSQUIDバイアスを設定する
2. 測定したい熱浴温度 $T_{\rm bath}$ に設定する
3. $I_{\rm bias}$ を2000µA流し、TES自身の発熱で完全に常伝導状態に転移させる
4. $I_{\rm bias}$ を1000µAまで徐々に小さくし、熱浴温度が安定するまで待つ
5. 熱浴温度が安定したら $I_{\rm bias}$ を徐々に下げながら出力電圧 $V_{\rm out}$ を測定する
6. $T_{\rm bath}$ を変えて2-5を同様に測定する

　$I_{\rm bias}$ と $V_{\rm out}$ の数値が測定できれば、式$(8),\ (9),\ (10)$から $I_{\rm TES},\ V_{\rm TES}$ が計算でき、I-V特性がわかります。その結果から、TESの諸パラメータの解析をしていきます。また、今回用いるデータにおける各定数のおおよその値は、それぞれ以下の通りです：

\begin{equation}
 \begin{aligned}
  M_{\rm in} &\sim 7.68\times10^{-11}\ {\rm H}\ , \\
  M_{\rm FB} &\sim 6.99\times10^{-11}\ {\rm H}\ , \\
  R_{\rm FB} &\sim 100\times10^3\ \Omega\ ,  \\
  R_{\rm sh} &\sim 3.9\times10^{-3}\ \Omega\ .
 \end{aligned}
\end{equation} \tag{17}

４ データ解析

　コードを見れば分かる方は、Google Colabで作成した以下の例を参照ください。この後に説明するコードも、以下のGoogle Colabに載せているものと同じです。

４.１ 解析準備

４.１.１ インストール

　解析には、非線形最小二乗法を用いたモデルフィッティングライブラリであるlmfitを使用するため、インストールしておきます。非線形最小二乗フィッティングは、測定データに最もよく適合する関数を見つける手法で、lmfitによってこれを簡単に行うことができます。具体的な使用法については、Ａ.１章にまとめています。
　pipはPythonのパッケージ管理ツールで、公開されているライブラリをダウンロードして、自分の環境にインストールするために使用します。行頭に!を付けることは、Google Colabのコードセルにおいてその行をPyhonコードではなくLinuxコマンドとして実行させるための合図になります。

!pip install lmfit

４.１.２ インポート

　使用するモジュールをインポートします。今回使用するモジュールは以下の通りです：

• glob：特定のパターンにマッチするファイルを取得するためのモジュール。
• numpy：配列などの数値計算を行うためのライブラリ。
• matplotlib.pyplot：グラフ描画用のモジュール。グラフ内のテキストや数式のフォントの指定もしておきます。
• matplotlib.cm：カラーマップを取得するためのモジュール。
• lmfit：scipy.optimizeを基盤とした、非線形最小二乗法によるフィッティングライブラリ。

　glob、numpy、matplotlibは、Google Colab環境に標準搭載されているため、追加のインストール作業は不要です。また、Pythonではasを使用することでライブラリやモジュールに短い別名を付けることができます。

import glob
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.family'] = 'serif'
plt.rcParams['mathtext.fontset'] = 'cm'
import matplotlib.cm as cm
import lmfit as lmf

４.１.３ データ取得

　GitHubで公開しているI-V測定データをGoogle Colab環境へダウンロードします。データは熱浴温度110mKから150mKまで5mK刻みで測定したI-Vデータです。ファイル内にはTESのバイアス電流値とそれに対するSQUIDの出力電圧値が並んでいます。
　wgetもpipと同様にLinuxコマンドで、指定したURLからファイルをネットワーク経由で取得します。このコードでは、リストに含まれる全ファイルをループ処理で順番にダウンロードし、Google Colabの作業スペース（カレントディレクトリ）に保存されます。

filename_list = [
  "IV_110mK_20230331.txt",
  "IV_115mK_20230331.txt",
  "IV_120mK_20230331.txt",
  "IV_125mK_20230331.txt",
  "IV_130mK_20230331.txt",
  "IV_135mK_20230331.txt",
  "IV_140mK_20230331.txt",
  "IV_145mK_20230331.txt",
  "IV_150mK_20230331.txt"
]

data_url = "https://raw.githubusercontent.com/Urchann/TES/main/TES_IV/IVdata_A5/"

for name in filename_list:
  !wget {data_url}{name}

４.１.４ 関数の定義

　解析に必要な関数の定義をします。関数は以下の通りです：

• TESに流れる電流の関数：式$(8)$を用います。
• TESの抵抗の関数：式$(9)$を用います。
• TESの両端にかかる電圧の関数：式$(10)$を用います。
• TESのJoule発熱の関数：式$(14)$を用います。
• TESの温度の関数：式$(16)$を用います。
• TESの熱伝導度の関数：式$(12)$を用います。
• TESの感度の関数：式$(1)$を用います。np.diffは、受け取った配列の隣接要素の差分を計算する関数です。 $[N_0,\ N_1,\ N_2,\ N_3]$ というデータがあれば、 $[N_1-N_0,\ N_2-N_1,\ N_3-N_2]$ を返します。このように、配列の長さは元のデータよりも一つ短くなるので、スライスで $T$ と $R$ の比の要素も揃える必要があります。スライスとは具体的に、データ $N=[N_0,\ N_1,\ N_2,\ N_3]$ に対して、N[1:]は $[N_1,\ N_2,\ N_3]$ 、N[:-1]は $[N_0,\ N_1,\ N_2]$ を抽出する処理になります。
• ファイル名から熱浴の温度を抽出する関数：splitは、文字列を分割するためのメソッドです。まず、split('/')[-1]でファイルパスのスラッシュ区切りで分割し、スライス[-1]、つまり末尾のファイル名を取り出します。次に、split('_')でファイル名をアンダースコア区切りで分解し、find('mK')で「mK」が含まれる要素のインデックスを特定します。最後に、その要素からスライス[:-2]、つまり二文字分の「mK」を取り除いてfloat型の数値に変換してTbathリストへ追加しています。

def Cal_Ites(Vout, Xi):
  """
  Calculate TES current.

  Parameters:
    Vout (list): Out put voltage of SQUID (unit: V).
    Xi  (float): Current-voltage conversion coefficient (unit: Ω).

  Returns:
    Ites (list): TES current (unit:A).
  """
  Ites = Vout / Xi
  return Ites


def Cal_Rtes(Itb, Ites, Rsh):
  """
  Calculate TES resistance.

  Parameters:
    Itb  (list): TES bias current (unit: A).
    Ites (list): TES current (unit: A).
    Rsh (float): Shunt resistance (unit: Ω).

  Returns:
    Rtes (list): TES resistance (unit: Ω).
  """
  Rtes = (Itb - Ites) * Rsh / Ites
  return Rtes


def Cal_Vtes(Rtes, Ites):
  """
  Calculate TES voltage.

  Parameters:
    Rtes (list): TES resistance (unit: Ω).
    Ites (list): TES current (unit: A).

  Returns:
    Vtes (list): TES voltage (unit: V).
  """
  Vtes = Rtes * Ites
  return Vtes


def Pfunc(Tb, Tc, G0, n):
  """
  Create a fitting function for TES heat generation.

  Parameters:
    Tbath (array): Heat bath temperature (unit: K).
    Tc    (float): Transition temperature (unit:K).
    G0    (float): Thermal conductivity as a heat bath temperature of 1mK (unit:W/K^n).
    n     (float): Temperature dependence of thermal conductivity (unit:None).

  Returns:
    pb    (array): TES heat generation (unit: W).
  """
  pb = (G0 / n)*(pow(Tc,n) - pow(Tbath,n))
  return pb


def Cal_Ttes(Tbath, G0, n, Ptes):
  '''
  Calculate TES temperature.

  Parameters:
    Tbath (array): Heat bath temperature (unit: K).
    G0    (float): Thermal conductivity as a heat bath temperature of 1mK (unit: W/K^n).
    n     (float): Temperature dependence of thermal conductivity (unit: None).
    Ptes  (array): TES heat generation (unit: W).

  Returns:
    Ttes  (array): TES temperature (unit: K).
  '''
  Ttes = ((Tbath**(n)) + ((n * Ptes) / G0))**(1 / n)
  return Ttes


def Cal_Gtes(G0, n, Ttes):
  '''
  Calculate TES thermal conductivity.

  Parameters:
    G0  (float): Thermal conductivity as a heat bath temperature of 1mK (unit: W/K^n).
    n   (float): Temperature dependence of thermal conductivity (unit: None).
    Ttes (list): TES temperature (unit: K).

  Returns:
    Gtes (list): TES thermal conductivity (unit: W/K).
  '''
  Gtes = G0 * (Ttes**(n-1))
  return Gtes


def Cal_alpha(Ttes, Rtes):
  """
  Calculate TES sensitivity.

  Parameters:
    Ttes   (list): TES temperature (unit: K).
    Rtes   (list): TES resistance (unit: Ω).

  Returns:
    alpha (float): TES sensitivity (unit: None).
  """
  dR = np.diff(Rtes)
  dT = np.diff(Ttes)
  dRdT = dR/dT
  T = (Ttes[1:] + Ttes[:-1])/2
  R = (Rtes[1:] + Rtes[:-1])/2
  alpha = dRdT * T / R
  return alpha


def makeTbath(filename_list):
  """
  Extract heat bath temperature values ​​from filenames.

  Parameters:
    filename_list (list): Measurement data files.

  Returns:
    Tbath         (list): Heat bath temperature (unit: mK).
  """
  Tbath = []
  for tb in filename_list:
    tb_ = tb.split('/')[-1]
    num = 0
    for i, name in enumerate(tb_.split('_')):
      if name.find('mK')!=-1:
        num = i
        break

    tb_key = tb_.split('_')[num]
    tbath  = float(tb_key[:-2])
    Tbath.append(tbath)

  Tbath = np.asarray(Tbath)
  return Tbath

４.１.５ 熱浴温度抽出

　makeTbath関数を用いて、各I-V測定データファイル名から測定した熱浴温度の数値を取得します。

Tbath = makeTbath(filename_list)

４.２ グラフ作成

４.２.１ データ計算

　データの取得ができたので、定義した関数を使って $I_{\rm TES},\ R_{\rm TES},\ V_{\rm TES}$ などの計算をします。まず、式$(17)$の通り、定数はここで与えておきます。さらに、各データの数値を格納するためのディクショナリを作成しておきます。
　それぞれの熱浴温度に対するデータをディクショナリに格納したいので、I-V測定データfilename_listを回します。このとき、４.１.５章で作ったTbathのリストから順に熱浴温度を引き出し、それを添え字として各ディクショナリに付ければ、それぞれの熱浴温度に対するデータ数値の格納が可能となります。各ディクショナリに入るデータは次の通りです：

• 読み取りデータdata：np.loadtxt関数によって、I-V測定データ（TXTファイル）に保存されている数値をNumPy配列として一括で読み込みます。また、データ点数もlen関数で取り出しています。
• 出力電圧 $V_{\rm out}$：出力電圧値はI-V測定データ（TXTファイル）の二列目に保存されているので、読み取ったデータからdata[:, 1]で抽出します。このとき、全体を最小出力電圧値で引いてオフセットをゼロにします。
• バイアス電流 $I_{\rm bias}$：バイアス電流値はI-V測定データ（TXTファイル）の一列目に保存されているので、読み取ったデータからdata[:, 0]で抽出します。このとき、測定されたバイアス電流の単位はµAなので、単位をAに変換するため全体を1e-6で割ります。
• TESに流れる電流 $I_{\rm TES}$：４.１.４章で定義したCal_Ites関数に $V_{\rm out}$ と $\Xi$ を代入します。
• TESの抵抗 $R_{\rm TES}$：４.１.４章で定義したCal_Rtes関数に $I_{\rm bias}$ と $R_{\rm sh}$ 、さらに計算した $I_{\rm TES}$ を代入します。
• TESの両端にかかる電圧 $V_{\rm TES}$：４.１.４章で定義したCal_Vtes関数に計算した $R_{\rm TES}$ と $I_{\rm TES}$ を代入します。
• TESのJoule発熱 $P_{\rm TES}$：式$(15)$の通り、TESのJoule発熱は $V_{\rm TES}I_{\rm TES}$ で書けるので、この計算をします。

　さて、TESのI-V特性グラフは３.３.２章で示した通り、横軸を $V_{\rm TES}$ 、縦軸を $I_{\rm TES}$ とします。しかし、計算に用いている $\Xi$ は $\mathcal L\gg1$ とした近似値であるため、超伝導領域（$T<T_{\rm c}$）において物理的にあり得ない負の傾きとなることがあります。この傾きを補正することで、物理的に妥当な特性曲線を得ます。
　まずは超伝導領域の傾き（${\rm d}I_{\rm TES}/{\rm d}V_{\rm TES}$）を算出し、その値をtiltとします。I-V測定データはバイアス電流の降順で記録されているため、データ配列の末尾付近が超伝導領域に対応します。そこで、データの末尾（最後から1番目）のバイアス電流は0Aであり除算エラーが生じるため除外し、最後の11番目から2番目までの計10点を用いて計算します。この計算には、４.１.４章の「TESの感度の関数」と同様にdiff関数を使います。
　続いて、算出したtiltが負であるかを判定して、負であればXiの値を1ずつ増加させながら $I_{\rm TES},\ R_{\rm TES},\ V_{\rm TES}$ を再計算し、tiltが正になるまでこの補正をループを繰り返します。
　最後に、filename_listの1回目のループで補正した $\Xi$ を用いて、各熱浴温度の $I_{\rm TES},\ R_{\rm TES},\ V_{\rm TES},\ P_{\rm TES}$ を一括して再度計算し、すべてのデータに対して補正を適用します。

Rsh = 3.9e-3
Min = 1.009e-10
Mfb = 8.912e-11
Rfb = 30.e3
Xi  = Rfb * Min / Mfb

data = {}
Vout = {}
Itb  = {}
Ites = {}
Rtes = {}
Vtes = {}
Ptes = {}

for i, filename in enumerate(filename_list):
  tb = Tbath[i]
  data     = np.loadtxt(filename)
  data_num = len(data[:, 0])
  Vout[tb] = data[:, 1] - data[:, 1].min()
  Itb[tb]  = data[:, 0] / 1e6
  Ites[tb] = Cal_Ites(Vout[tb], Xi)
  Rtes[tb] = Cal_Rtes(Itb[tb], Ites[tb], Rsh)
  Vtes[tb] = Cal_Vtes(Rtes[tb], Ites[tb])
  Ptes[tb] = Ites[tb] * Vtes[tb]

  if i == 0:
    tilt = np.average(np.diff(Ites[tb][data_num-11:data_num-1]))/np.average(np.diff(Vtes[tb][data_num-11:data_num-1]))
    print(f'T    = {tb} mK')
    print(f'Xi   = {Xi} Ohm')
    print(f'tilt = {tilt}')
    j = 0
    while tilt < 0:
      j  += 1
      Xi += 1
      Ites[tb] = Cal_Ites(Vout[tb], Xi)
      Rtes[tb] = Cal_Rtes(Itb[tb], Ites[tb], Rsh)
      Vtes[tb] = Cal_Vtes(Rtes[tb], Ites[tb])
      tilt     = np.average(np.diff(Ites[tb][data_num-11:data_num-1]))/np.average(np.diff(Vtes[tb][data_num-11:data_num-1]))

    print(f'{j} modified Xi   = {Xi} Ohm')
    print(f'{j} modified tilt = {tilt}')

  Ites[tb] = Cal_Ites(Vout[tb], Xi)
  Rtes[tb] = Cal_Rtes(Itb[tb], Ites[tb], Rsh)
  Vtes[tb] = Cal_Vtes(Rtes[tb], Ites[tb])
  Ptes[tb] = Ites[tb] * Vtes[tb]

４.２.２ I-V特性グラフ

　まずはI-V特性グラフを見てみます。forループで各熱浴温度に対する $I_{\rm bias}$ と $V_{\rm out}$ の関係、および $V_{\rm TES}$ と $I_{\rm TES}$ の関係を一つのグラフにそれぞれ重ねています。プログラムの詳解については、Matplotlibの知識しか使っていないため本記事では割愛します。
　左の図が横軸 $I_{\rm bias}$ 、縦軸 $V_{\rm out}$ の素データで、右の図が横軸 $V_{\rm TES}$ 、縦軸 $I_{\rm TES}$ のデータの特性グラフになっています。

_, ax = plt.subplots(1, 2, figsize=(15, 5))
for i in Tbath:
  ax[0].plot(Itb[i]*1e6, Vout[i], '.', label=r'$T_{\rm bath}=$'+str(int(i))+r'$\ {\rm mK}$')
  ax[1].plot(Vtes[i]*1e6, Ites[i]*1e6, '.', label=r'$T_{\rm bath}=$'+str(int(i))+r'$\ {\rm mK}$')

ax[0].set_title(r'$V_{\rm out}\ {\rm vs.}\ I_{\rm TES}$', fontsize=16)
ax[0].set_xlabel(r'$I_{\rm TES}\ ({\rm\mu A})$', fontsize=14)
ax[0].set_ylabel(r'$V_{\rm out}\ ({\rm V})$', fontsize=14)
ax[0].legend(fontsize=12, frameon=False, ncols=2)
ax[0].grid(ls=':', alpha=.5)
ax[1].set_title(r'$I_{\rm TES}\ {\rm vs.}\ V_{\rm TES}$', fontsize=16)
ax[1].set_xlabel(r'$V_{\rm TES}\ ({\rm V})$', fontsize=14)
ax[1].set_ylabel(r'$I_{\rm TES}\ ({\rm\mu A})$', fontsize=14)
ax[1].legend(fontsize=12, frameon=False, ncols=2)
ax[1].grid(ls=':', alpha=.5)
plt.show()

４.２.３ P-R特性グラフ

　次にP-R特性グラフも見てみましょう。print文でRtesの数値を確認するとわかりますが、最後にnanの文字があると思います。４.２.１章のコードを実行するとRuntimeWarning: invalid value encountered in divideのような表示があるはずで、これは $I_{\rm TES}$ の中にあるゼロで割っている要素があるためです。今回のプログラムの実行にはさほど影響しませんが、扱いには気を付ける必要があります。
　具体的には、for文内のrtes = Rtes[tb][Rtes[tb]==Rtes[tb]]のスライス部分で、NaN（Not a Number）以外の要素の呼び出しを行います。ここでは、NaN値の どの値とも等しくない（自身とも等しくない） という性質を用いて、[リスト=リスト]という一見無意味に思えるスライスによって、NaN値のみが除外されたデータに焼き直しているのです。
　このようにして得た $R_{\rm TES}$ の配列から最小値と最大値を取り出し、 $R_{\rm TES}$ を正規化した値を横軸、TESJouleの発熱 $P_{\rm TES}$ を縦軸にしてグラフにしています。その他のプログラムは先ほどと同様にMatplotlibの知識しか使っていないため、省略します。

_, ax = plt.subplots(1, 1, figsize=(6, 4))
for tb in Tbath:
  rtes    = Rtes[tb][Rtes[tb]==Rtes[tb]]
  rtesmax = rtes.max()
  rtesmin = rtes.min()
  Rn      = rtesmax - rtesmin
  ax.plot((Rtes[tb]-rtesmin)/Rn, Ptes[tb]*1e9, '.', label=r'$T_{\rm bath}=$'+str(int(tb))+r'$\ {\rm mK}$')

ax.set_title(r'$P_{\rm TES}\ {\rm vs.}\ R_{\rm TES}$', fontsize=16)
ax.set_xlabel(r'$R_{\rm TES}/R_{N}\ ({\rm\mu A})$', fontsize=14)
ax.set_ylabel(r'$P_{\rm TES}\ ({\rm nW})$', fontsize=14)
ax.legend(fontsize=12, frameon=False, ncols=2)
ax.grid(ls='--', alpha=.5)
ax.set_ylim(0., 0.1)
# plt.xlim(0., 300)
plt.show()

４.２.４ フィッティング

　さて、ここでTESのパラメータを同定します。式$(15)$の通り、TESのJoule発熱と熱浴温度の間には

P_{\rm TES} = \frac{G_0}{n}(T_{\rm TES}^n-T_{\rm bath}^n)

の関係がありました。TESの温度を転移温度 $T_{\rm c}$ として、測定データに対してこの関数でフィッティングを行うことで、 $G_0,\ n,\ T_{\rm c}$ の各パラメータを算出します。ただし、 $y=ax^n-b$ という冪函数のフィットは、求めたい変数（$G_0,\ n,\ T_{\rm c}$）と自由度がともに3つで原理的には決まりそうに見えるものの、関数の範囲が狭い（直線に近い）ので $G_0$ と $n$ が縮退しやすくなります。そのため、データによっては $n$ を3または4などに固定しても良いです。
　３.３.２章で述べたように、転移温度付近では $P_{\rm TES}$ が一定値になります。４.２.３章の結果を見れば、該当するのは $R_{\rm TES}/R_N=0.1\sim0.7$ の辺りであることがわかります。なるべく $P_{\rm TES}$ が安定している範囲で考えたいので、今回は適当に $R_{\rm TES}/R_N=0.2\sim0.3$ を採用します。
　まず、for文内で各熱浴温度について途中までは４.２.３章と同様の計算をし、 $P_{\rm TES}$ の値として $R_{\rm TES}/R_N=0.2\sim0.3$ の数値を取り出します。具体的には、(0.2<rtes) & (rtes<0.3)とすることでrtesの値が0.2より大きく0.3より小さい部分はTrue、それ以外はFalseと変換されるため、このフィルタによってスライス指定しています。
　続いて、取得した転移温度付近の $P_{\rm TES}$ の中で平均をとり、さらに計算した際のデータの不偏標準偏差から誤差を計算します。不偏標準偏差の計算にはNumPyのstd関数を使います。蛇足ですが、Tb=Tbath/1e3は熱浴温度の単位をmKからKに変換しているだけです。
　転移温度付近におけるフィッティング対象のデータPbとその誤差Pberrの準備ができたので、さっそくフィッティングを行っていきましょう。以下にlmfitモジュールから行っているプログラムを少しまとめます：

• model = lmf.Model()：フィッティングモデル関数の定義です。今回は４.１.４章で定義したPfunc関数を代入して使用します。
• model.make_params()：初期パラメータの生成を呼び出し、詳細な情報が出力されるように設定します。verbose=Trueは、今回のプログラムを実行すれば次のような出力が得られるだけで、無くても問題ありません。

  - Adding parameter "Tc"
  - Adding parameter "G0"
  - Adding parameter "n"
  

• model.set_param_hint()：フィッティングに用いるパラメータの初期値、範囲、拘束条件などを事前に設定します。今回のフィッティングにおけるパラメータは $G_0,\ n,\ T_{\rm c}$ で、それぞれの考えられる最小値と最大値の設定をここで行っています。
• result = model.fit()：フィッティングの実行。引数には一番最初にフィッティング対象のデータ、続いてモデル関数の引数とその初期パラメータ値を順に入れます。フィッティングでは、初期パラメータ値を中心にmodel.set_param_hint()で指定した拘束条件の中で動かし、最適パラメータを探しています。
• print(result.fit_report())：フィッティング結果を解析し、表示させます。今回のプログラムを実行すれば次のような出力が得られます。

  [[Model]]
      Model(Pfunc)
  [[Fit Statistics]]
      # fitting method   = leastsq
      # function evals   = 152
      # data points      = 9
      # variables        = 3
      chi-square         = 1.2775e-25
      reduced chi-square = 2.1291e-26
      Akaike info crit   = -529.652725
      Bayesian info crit = -529.061051
      R-squared          = 1.00000000
  [[Variables]]
      Tc:  0.15526184 +/- 2.8240e-04 (0.18%) (init = 0.29)
      G0:  3.0544e-08 +/- 8.5904e-09 (28.12%) (init = 1e-08)
      n:   3.00000640 +/- 0.17114333 (5.70%) (init = 3.5)
  [[Correlations]] (unreported correlations are < 0.100)
      C(G0, n)  = +0.9998
      C(Tc, G0) = -0.8297
      C(Tc, n)  = -0.8189
  

• result.best_values[]：[]内にフィッティングパラメータを代入すると、それぞれの最適フィット値が返されます

　最後に、chi2_minでフィッティング結果について定量的な精度評価をカイ二乗検定により行います。その後のプロットに関するMatplotlibのプログラムについては、これまで同様省略します。

Pb    = []
Pberr = []
for tb in Tbath:
  rtes_   = Rtes[tb][Rtes[tb]==Rtes[tb]]
  rtesmax = rtes_.max()
  rtesmin = rtes_.min()
  Rn      = rtesmax - rtesmin
  rtes    = Rtes[tb]/Rn

  filter_Ptes = (.2<rtes) & (rtes<.3)

  pb = np.average(Ptes[tb][filter_Ptes])
  Pb.append(pb)

  perr = np.std(Ptes[tb][filter_Ptes], ddof=1)
  Pberr.append(perr)

Tb    = Tbath/1e3
Pb    = np.asarray(Pb)
Pberr = np.asarray(Pberr)

model = lmf.Model(Pfunc)
model.make_params(verbose=True)
model.set_param_hint('Tc',min=.1, max=.32)
model.set_param_hint('G0',min=1.e-15, max=1.e2)
model.set_param_hint('n', min=3, max=4)
result = model.fit(Pb, Tbath=Tb, Tc=.29, G0=1.e-8, n=3.5)
print(result.fit_report())
Tc = result.best_values['Tc']
G0 = result.best_values['G0']
n  = result.best_values['n']

chi2_min = sum(((Pb - Pfunc(Tbath=Tb, Tc=Tc, G0=G0, n=n))/Pberr)**2)
print(f"χ² = {chi2_min}")

_, ax = plt.subplots(1, 1, figsize=(6, 4))
ax.errorbar(Tb*1e3, Pb*1e12, yerr = Pberr*1e12, fmt='ok', capsize=4, ms=5)
tb_min = Tbath.min()/1e3 - 0.010
tb_max = Tbath.max()/1e3 + 0.010
t = np.linspace(tb_min, tb_max, 100)
ax.plot(t*1e3, Pfunc(t, Tc, G0, n)*1e12, 'r--')
ax.set_title(r'$P_{\rm bath}\ {\rm vs.}\ T_{\rm bath}$', fontsize=16)
ax.set_xlabel(r'$T_{\rm bath}\ ({\rm mK})$', fontsize=14)
ax.set_ylabel(r'$P_{\rm bath}\ ({\rm pW})$', fontsize=14)
ax.legend(['data',rf'$T_{{\rm c}}={Tc:.3}\ {{\rm K}},\ G_0={G0*1e9:.3}\ {{\rm nW/K}},\ n={n:.2f}$'], fontsize=12, frameon=False)
ax.grid(ls='--', alpha=.5)
plt.show()

４.２.５ R-T特性グラフ

　$G_0$ と $n$ の最適フィット値が求められたので、４.１.４章で定義したCal_Ttes関数により $T_{\rm TES}$ を計算することができます。計算した $T_{\rm TES}$ と $R_{\rm TES}$ の特性を見れば、３.１ 温度計の感度で見たグラフと類似した結果が得られます。

#Calculation T_TES
for i in Tbath:
  Ttes[i] = Cal_Ttes(i/1e3, Go, n, Ptes[i])
  plt.plot(Ttes[i]*1e3, Rtes[i], '.',label='Tbath='+str(int(i))+'mK')

plt.title(r'$\rm R_{TES}\ vs.\ T_{TES}$')
plt.xlabel(r'$\rm T_{TES}\ (mK)$', fontsize=14)
plt.ylabel(r'$\rm R_{TES}\ (\Omega)$', fontsize=14)
plt.legend(fontsize=10,frameon=False, ncol=1,loc=2)
plt.grid(ls='--')

４.２.６ TESの熱伝導度

　$G$ についても計算してみましょう。４.１.４章で定義したCal_Gtes関数に各パラメータを代入して、 $T_{\rm TES}$ との関係を見ます。式$(12)$より、 $G$ と $T_{\rm TES}$ の関係は

G = G_0T_{\rm TES}^{n-1}

で、 $n$ の最適フィット値は3.00でだったので、グラフは二次関数となります。

Gtes = {}

_, ax = plt.subplots(1, 1, figsize=(6, 4))
for tb in Tbath:
  Gtes[tb] = Cal_Gtes(G0, n, Ttes[tb])
  ax.plot(Ttes[tb]*1e3, Gtes[tb], '.',label='Tbath='+str(int(i))+'mK')

ax.set_title(r'$G_{\rm TES}\ {\rm vs.}\ T_{\rm TES}$', fontsize=16)
ax.set_xlabel(r'$T_{\rm TES}\ ({\rm mK})$', fontsize=14)
ax.set_ylabel(r'$G_{\rm TES}\ ({\rm W/K})$', fontsize=14)
ax.legend(fontsize=12, frameon=False, ncol=1, loc=2)
ax.grid(ls='--', alpha=.5)
plt.show()

４.２.７ TESの感度

　$\alpha$ についても計算します。４.１.４章で定義したCal_alpha関数に、ここでは熱浴温度120mKにおけるパラメータを代入して、バイアス電流$I_{\rm bias}$との関係を見ます。Cal_alpha関数での定義で $T$ と $R$ を隣接する要素の平均値として取り出したため、比較する $I_{\rm bias}$ の値も隣接する要素の平均値を取り出しています。
　比較する量が$I_{\rm bias}$であるのは、実験において自由に可変な量がバイアス電流のため、このようなグラフを作ることによって、バイアス電流がどの値で超伝導転移端と対応するのかを見たいという理由があります。論文で結果を書く際には、横軸を $R_{\rm TES}/R_{\rm N}$ にします。
　定義から、$\alpha$ の値は超伝導遷移端で大きくなるのですが、R-T特性グラフの転移温度付近のデータ点をよく見てみると、隣接するデータ点が完全な増加傾向になっておらず、所々水平に近い部分があります。そのため、下のグラフにおいても $\alpha$ の値が不規則にピークを持った形をしています。

alpha = Cal_alpha(Ttes[120], Rtes[120])
x = (Itb[120][1:] + Itb[120][:-1])/2

plt.plot(x*1e6, alpha)
plt.title(r'$\rm alpha$')
plt.xlabel(r'$\rm I_{bias}\ (\mu\rm A)$', fontsize=14)
plt.ylabel(r'$\alpha$', fontsize=14)
plt.ylim(0, 200)
plt.xlim(150,300)

４.３ コントアの計算

　ここまで $G_0,\ n,\ T_{\rm c}$ の最適フィット値を使ってきましたが、フィッティングによって得られた数値は、分散が最小になるよう計算して収束したときの結果に過ぎず、統計上重要なのはその信頼範囲です。そこで、 $G_0,\ n,\ T_{\rm c}$ について二つずつ各組み合わせで信頼区間（コントア）を計算してみます。
　ちなみに、コントア（contour）とは「輪郭」「等高線」といった意味があります。

４.３.１ 関数の定義

　まずはコントアの計算に必要な関数を導入します：

• 信頼範囲の定義：今回、コントア計算関数の自由度はフィッティング変数のうち二つであるため、自由度2のカイ二乗分布における累積確率の値からパーセント点を求めます。累積確率は68.27%、90%、99%の三点のP値を考え、ここのサイトから取ってきた値を使います。「累積モード」は下側累積$\rm P$に、「自由度」は2にセットして、「累積確率」は0.6827、0.9、0.99のそれぞれを入れれば取得できます。

• コントア関数の定義：名前通りですが、CalCont_G0vsn関数は $G_0$ と $n$ のコントアを計算する関数、CalCont_Tcvsn関数は $T_{\rm c}$ と $n$ のコントアを計算する関数、CalCont_TcvsG0関数は $T_{\rm c}$ と $G_0$ のコントアを計算する関数です。それぞれの関数ではほとんど同じ計算をするので、以降はCalCont_G0vsn関数についてのみ解説します
  　関数の引数は、熱浴温度Tb、各熱浴温度における転移温度付近の $P_{\rm TES}$ Pbとその誤差Perr、 $T_{\rm c}$ の最適フィット値Tcbest、カイ二乗値chi2_min、 $n$ の計算区間配列xn、 $G_0$ の計算区間配列ygです。また、結果を保存するときのファイル名は設定せず、name=Noneとしておきます。

• メッシュグリッドの作成：CalCont_G0vsn関数では、 $n$ と $G_0$ の格子配列を作成します。この意義は、XY平面全域にわたる信頼度計算のための格子点生成にあります。np.meshgrid関数の操作は下のようなイメージです（下図はここのサイトから引用）。

  

  　二次元配列のXX, YYおよびYY.T（YYの転置をとり、列成分の要素を配列にする）は以下のように生成されることを抑えておきましょう：

  import numpy as np
  
  x = np.arange(0, 6, 1)
  y = np.arange(0, 4, 1)
  XX, YY = np.meshgrid(x, y)
  
  XX, YY, YY.T
  '''
  (array([[0, 1, 2, 3, 4, 5],
          [0, 1, 2, 3, 4, 5],
          [0, 1, 2, 3, 4, 5],
          [0, 1, 2, 3, 4, 5]]),
   array([[0, 0, 0, 0, 0, 0],
          [1, 1, 1, 1, 1, 1],
          [2, 2, 2, 2, 2, 2],
          [3, 3, 3, 3, 3, 3]]))
   array([[0, 1, 2, 3],
          [0, 1, 2, 3],
          [0, 1, 2, 3],
          [0, 1, 2, 3],
          [0, 1, 2, 3],
          [0, 1, 2, 3]]))
  '''
  

• カイ二乗値の計算：まず、計算したカイ二乗値を格納するための空配列を用意しておき、その名前はchi2_Fix_nとします。また、あとでCalCont_G0vsn関数を使う際に、引数のg0をnW/K^nの単位で代入するため、yg配列の要素g0の単位を$\rm W/K$に直しています。次に $T_{\rm c}$ を最適フィット値で固定して、xnとygの中で $n,\ G_0$ を動かし、カイ二乗値を計算します。for文が連なっており少し複雑ですが、

  ({\rm xn}_0,\ {\rm yg}_0),\ ({\rm xn}_0,\ {\rm yg}_1),\ \cdots,\ ({\rm xn}_0,\ {\rm yg}_n),\ ({\rm xn}_1,\ {\rm yg}_0),\ ({\rm xn}_1,\ {\rm yg}_1),\ \cdots({\rm xn}_m,\ {\rm yg}_n)
  

  の組み合わせ順で計算されています。蛇足ですが、添え字が0から始まっているのは、enumerate関数ではインデックス番号がデフォルトで0から開始されるためです。
  　print文が少々複雑ですが、これはループの進行状況をリアルタイムで標準出力する（1000回ループするような処理があったとき、進行状況を「$1/1000$」→「$2/1000$」のように表示させる）プログラムで、無くても問題ありません。ただ、せっかくなので簡単に説明しておきます：

  • f文字列
    　print文に入れる文字列の先頭にfまたはFを付け、その中で変数や式を中カッコで囲むことで、その値を文字列に埋め込むことができます。以下の具体例を見た方が早いでしょう：

    name = "Taro"
    age  = 20
    
    print('{name} is {age} years old.')
    # {name} is {age} years old.
    
    print(f'{name} is {age} years old.')
    # Taro is 20 years old.
    

  • \r（キャリッジリターン）
    　次に続く文字が現在の行を上書きするという作用を持ちます。以下で具体的に見てみましょう（print('123456789\rABC')を実行してみると、キャリッジリターンの特性がさらにわかると思います）：

    print('Hello\nWorld')
    #Hello
    #World
    
    print('Hello\rWorld')
    #World
    

    　ここで注意したいのは、テキスト内で\rを使うと、カーソルが同じ行の先頭に戻り、直後の文字列が上書きされるのであって、今回のプログラムようにfor文でループさせた文字列では自動的に改行が入るため、文字列は上書きされず、\rを入れたときと入れないときとでプログラムの実行結果は変化しません（試しに、\rおよびをend=''を消した場合とend=''のみを消した場合でプログラムを実行してみてください）。つまり、キャリッジリターンだけでは、進行状況をリアルタイムで標準出力するには不十分で、次に紹介するendオプションが役に立ちます
  • end=''（endオプション）
    　上で述べたように、print文では行末にデフォルトで改行が入ってしまいます。これを防ぐために、endオプションにより行末の表示文字列をデフォルトの改行文字\nから何もなし''に指定することができます。以下で具体的に見てみましょう：

    print('abc')
    print('def')
    #abc
    #def
    
    print('abc', end='')
    print('def', end='')
    #abcdef
    

    　これを使えば、for文でループさせた文字列の改行を消せて、\r文字列\r文字列\r文字列のようになるため、キャリッジリターンとの組み合わせで進行状況をリアルタイムで標準出力することができます。
  • flush=True
    　print文では、効率化の観点から出力がバッファリングされることがあります。せっかく進行状況をリアルタイムで確認できるプログラムを組んだというのにバッファリングされては意味がありません。そこで、print関数の引数にflush=Trueと書くことにより、即座に出力を画面に表示させることを強制できます。ただ、今回のプログラムに限っては、このオプションは不要かもしれません。
    　以上を踏まえ、改めてprint(f'\r{i+1}/{len(Xn[0])},{j+1}/{len(YG0.T[0])}',end='',flush=True)を見てみましょう。すでに述べた通り、enumerate関数ではインデックス番号が0から開始されるため、{i+1}, {j+1}として進行状況が1から開始されるようにしています。また、進行状況の分母と分子はxnとygの要素数を入れれば良く、メッシュグリッドの作成における図で見たことを思い出せば、配列xnの要素数はlen(Xn[0])で取り出せ、配列ygの要素数は、YG0の転置（YG0.T）をとれば、その一つ目の配列の長さlen(YG0.T[0])で取り出せることがわかると思います。

• カイ二乗値の計算結果の格納：さて、計算したカイ二乗値の結果は

  \begin{pmatrix}
  {\rm chi}2({\rm xn}_0,\ {\rm yg}_0)
  &
  {\rm chi}2({\rm xn}_0,\ {\rm yg}_1)
  &
  \cdots
  &
  {\rm chi}2({\rm xn}_0,\ {\rm yg}_n)
  \\
  {\rm chi}2({\rm xn}_1,\ {\rm yg}_0)
  &
  {\rm chi}2({\rm xn}_1,\ {\rm yg}_1)
  &
  \cdots
  &
  {\rm chi}2({\rm xn}_1,\ {\rm yg}_n)
  \\
  \vdots & \vdots & \ddots & \vdots
  \\
  {\rm chi}2({\rm xn}_m,\ {\rm yg}_0)
  &
  {\rm chi}2({\rm xn}_m,\ {\rm yg}_1)
  &
  \cdots
  &
  {\rm chi}2({\rm xn}_m,\ {\rm yg}_n)
  \end{pmatrix}
  

  のように格子配列にして格納したいです。そこで使うのが、指定した配列同士を横に重ねるnp.hstack関数と、指定した配列同士を縦に重ねるnp.vstack関数です。for j, g0 in numerate(yg)の中では、上に示した格子配列の行成分を作ることができるため、先に生成しておいた空配列chi2_Fix_nとカイ二乗値の計算結果chi_contをnp.hstack関数で横に重ねます。またfor i, nn in numerate(xn)の中では、enumerate関数のインデックス番号が1以降のとき、横に重なった配列をnp.vstack関数で縦に重ねて、上のような格子配列を生成することができます。

• ゼロ点調整：カイ二乗値の計算結果chi_contの最小値をゼロにするため、４.２.４章で計算した最小カイ二乗値chi_minにより全体を引いておきます。

• コントア図の描画：ax.contour()によりコントア図を作ることができます。引数には通常、X座標を表す一次元配列 $(M,\ 1)$ 、Y座標を表す一次元配列 $(1,\ N)$ 、各座標に対する二次元配列 $(M,\ N)$ 、等高線の間隔levelsを指定します。levelsには、各座標に対する二次元配列 $(M,\ N)$ の値が $a,\ 2a,\ 3a,\cdots$ となるコントアを作成したいのであればlevels=aとすればよく、任意の値で指定したいのであればlevels=[a,b,c,...]とすれば良いです。
  　CalCont_G0vsn関数では、X座標は $n$ にするのでXn、Y座標は $G_0$ にするのでYGoです。したがって、各座標に対する二次元配列はchi2_arrayの転置をとったものになります。また、等高線の間隔は信頼範囲cslevelsです。さらに、clabel(fmt='%1.1f', fontsize=14)によって等高線にフォントサイズ14、一桁の小数点以下を持つ浮動小数点数のラベルを付けます。

• フィッティング値の描画：CalCont_Govsn関数では、フィッティング値$n,\ G_0$のプロットもしてきます。マーカーはサイズ$10$の星形状にします（plt.plot(n, Go*1e9, 'r*', ms=10)）。

• データの保存：関数の引数でnameが存在する場合、たとえばnp.savetxt(f'G0vsn_X_{name}.txt',X)ならnameに定義した名前が入り、G0vsn_X_{name}.txtというファイル名で、開いているマイドライブにXnデータが保存されます。今回はname=Noneと定義したので、デフォルトではデータが保存されません。

CRange68 = 2.295815160785974337606
CRange90 = 4.605170185988091368036
CRange99 = 9.210340371976182736072
cslevels = [CRange68,CRange90,CRange99]


def CalCont_G0vsn(Tb, Pb, Perr, Tcbest, chi2_min, xn, yg, name=None):
  Xn, YG0 = np.meshgrid(xn, yg)
  for i, nn in enumerate(xn):
      chi2_Fix_n = np.array([])
      for j, g0 in enumerate(yg):
          g0 = g0*1e-9
          chi2_cont = sum(((Pb - Pfunc(Tbath=Tb, Tc=Tc, G0=g0, n=nn))/Perr)**2)
          print(f'\rCalculation contour: G0 vs n  --> {i+1}/{len(Xn[0])}, {j+1}/{len(YG0.T[0])}', end='', flush=True)
          # print(f'\r{e+1}/{len(Xn[0])}',end='',flush=True)

          chi2_Fix_n = np.hstack((chi2_Fix_n, chi2_cont))

      if i == 0:
        chi2_array = chi2_Fix_n
      else:
        chi2_array = np.vstack((chi2_array, chi2_Fix_n))

  print()
  chi2_array -= chi2_min

  _, ax = plt.subplots(1, 1, figsize=(6, 4))
  cont = plt.contour(Xn, YG0, chi2_array.T, levels=cslevels)
  cont.clabel(fmt='%1.1f', fontsize=14)
  ax.set_title(r'$G_0\ {\rm vs.}\ n$'+' Contour Plot', fontsize=16)
  ax.plot(n, G0*1e9, 'r*', ms=10)
  ax.set_xlabel(r'$n$', fontsize=14)
  ax.set_ylabel(r'$G_0\ ({\rm nW/K^n})$', fontsize=14)
  if name:
      np.savetxt(f'G0vsn_X_{name}.txt', Xn)
      np.savetxt(f'G0vsn_Y_{name}.txt', YG0)
      np.savetxt(f'G0vsn_Z_{name}.txt', chi2_array)
      # np.savetxt(f'BestFit_{name}.txt', bestfit_para)


def CalCont_Tcvsn(Tb, Pb, Perr, Gobest, chi2_min, xn, ytc, name=None):
  Xn, YTc = np.meshgrid(xn, ytc)
  for i, nn in enumerate(xn):
      chi2_Fix_n = np.array([])
      for j, tc in enumerate(ytc):
          tc = tc*1e-3
          chi2_cont = sum(((Pb - Pfunc(Tbath=Tb, Tc=tc, G0=G0, n=nn))/Perr)**2)
          print(f'\rCalculation contour: Tc vs n  --> {i+1}/{len(Xn[0])},{j+1}/{len(YTc.T[0])}',end='',flush=True)
          # print(f'\r{e+1}/{len(X[0])}',end='',flush=True)

          chi2_Fix_n = np.hstack((chi2_Fix_n, chi2_cont))

      if i == 0:
        chi2_array = chi2_Fix_n
      else:
        chi2_array = np.vstack((chi2_array, chi2_Fix_n))

  print()
  chi2_array -= chi2_min

  _, ax = plt.subplots(1, 1, figsize=(6, 4))
  cont = ax.contour(Xn, YTc, chi2_array.T, levels=cslevels)
  cont.clabel(fmt='%1.1f', fontsize=14)
  ax.set_title(r'$T_{\rm c}\ {\rm vs.}\ n$'+' Contour Plot', fontsize=16)
  ax.plot(n, Tc*1e3, 'r*', ms=10)
  ax.set_xlabel(r'$n$', fontsize=14)
  ax.set_ylabel(r'$T_{\rm c}\ ({\rm mK})$', fontsize=14)
  if name:
      np.savetxt(f'Tcvsn_X_{name}.txt', Xn)
      np.savetxt(f'Tcvsn_Y_{name}.txt', YTc)
      np.savetxt(f'Tcvsn_Z_{name}.txt', chi2_array)
      # np.savetxt(f'BestFit_{name}.txt', bestfit_para)


def CalCont_TcvsG0(Tb, Pb, Perr, nbest, chi2_min, xg, ytc, name=None):
  XG0, YTc = np.meshgrid(xg, ytc)
  for i, g0 in enumerate(xg):
      chi2_Fix_n = np.array([])
      g0 = g0*1e-9
      for j, tc in enumerate(ytc):
          tc = tc*1e-3
          chi2_cont = sum(((Pb - Pfunc(Tbath=Tb, Tc=tc, G0=g0, n=n))/Perr)**2)
          print(f'\rCalculation contour: Tc vs G0 --> {i+1}/{len(XG0[0])},{j+1}/{len(YTc.T[0])}',end='',flush=True)
          # print(f'\r{e+1}/{len(X[0])}',end='',flush=True)

          chi2_Fix_n = np.hstack((chi2_Fix_n, chi2_cont))

      if i == 0:
        chi2_array = chi2_Fix_n
      else:
        chi2_array = np.vstack((chi2_array, chi2_Fix_n))

  print()
  chi2_array -= chi2_min

  _, ax = plt.subplots(1, 1, figsize=(6, 4))
  cont = ax.contour(XG0, YTc, chi2_array.T, levels=cslevels)
  cont.clabel(fmt='%1.1f', fontsize=14)
  ax.set_title(r'$T_{\rm c}\ {\rm vs.}\ G_0$'+' Contour Plot', fontsize=16)
  ax.plot(G0*1e9, Tc*1e3, 'r*', ms=10)
  ax.set_xlabel(r'$G_0\ ({\rm nW/K^n})$', fontsize=14)
  ax.set_ylabel(r'$T_{\rm c}\ ({\rm mK})$', fontsize=14)
  if name:
      np.savetxt(f'TcvsGo_X_{name}.txt', XG0)
      np.savetxt(f'TcvsGo_Y_{name}.txt', YTc)
      np.savetxt(f'TcvsGo_Z_{name}.txt', chi2_array)
      # np.savetxt(f'BestFit_{name}.txt', bestfit_para)

４.３.２ 計算の実行

　４.３.１章で与えたCalCont_G0vsn関数、CalCont_Tcvsn関数、およびCalCont_TcvsG0関数のそれぞれの引数に変数を代入します。 $G_0,\ n,\ T_{\rm c}$ のそれぞれのメッシュ配列範囲は任意に変えましょう。

CalCont_G0vsn(Tb, Pb, Pberr, Tc, chi2_min, xn=np.arange(2.5, 3.5, 0.01), yg=np.arange(0, 60, 0.1))
CalCont_Tcvsn(Tb, Pb, Pberr, G0, chi2_min, xn=np.arange(2.95, 3.05, 0.01), ytc=np.arange(154, 156.5, 0.01))
CalCont_TcvsG0(Tb, Pb, Pberr, n, chi2_min, xg=np.arange(20, 40, 1), ytc=np.arange(153, 157, 0.1))

５ まとめ

　本記事では、TESのI-V測定で得られるバイアス電流とSQUID出力電圧から、TESに流れる電流、TESの両端にかかる電圧、TESの抵抗値、TESのJoule発熱を求める一連の解析方法を紹介しました。
　解析の流れを整理すると、次のようになります。

1. SQUIDの出力電圧 $V_{\rm out}$ をTESに流れる電流 $I_{\rm TES}$ に変換する。
2. シャント抵抗を含むバイアス回路の関係式から、TESの両端にかかる電圧 $V_{\rm TES}$ とTESの抵抗 $R_{\rm TES}$ を求める。
3. I-V特性およびP-R特性から、超伝導領域・転移領域・常伝導領域を確認する。
4. 転移領域のTESのJoule発熱 $P_{\rm TES}$ と熱浴温度 $T_{\rm bath}$ の関係をフィッティングし、転移温度 $T_{\rm c}$ 、熱輸送の係数 $G_0$ 、温度依存性の指数 $n$ を推定する（ただし、$G_0$ と $n$ は縮退しやすい）。
5. 得られたパラメータから、R-T特性、TESの熱伝導度 $G$ 、TESの感度 $\alpha$ を評価する。
6. パラメータ間のコントアを計算し、最適フィット値だけでなく、その不確かさや相関も確認する。

　このように、I-V測定は単に電流と電圧の関係を描くための測定ではありません。測定された電気信号を、回路モデルと熱モデルを介して解釈することで、TESの転移温度、熱浴との結合、温度に対する感度といった、デバイスの動作を支配する物理量を引き出すことができます。
　一方で、得られる結果は電流電圧変換係数 $\Xi$ 、シャント抵抗 $R_{\rm sh}$ 、解析に用いる転移領域の選び方、測定点のばらつきなどにも依存します。したがって、フィッティングによる最適値だけを見るのではなく、元のI-V曲線との整合性、パラメータ間の相関、信頼区間まで含めて評価することが重要です。

補遺Ａ

　本記事のコードを実装する上でポイントとなった手法について、ここで補足します。

Ａ.１ lmfitの使い方

　lmfitにより、フィッティングモデルの定義、パラメータの設定、フィッティングの実行、および結果の解析を簡単に行うことができます。端的には、あるパラメータを持ったフィッティングモデル関数を用意して、解析に使うデータとうまく合うような最適パラメータを見つけるものです。以下に基本的な解析手順を説明します：

1. lmfitのインストール：まずは、pipを使ってlmfitをインストールしておきましょう。
2. データの準備：フィッティングするデータを用意します。
3. モデルの定義：フィッティングに使用するモデル関数を定義します。モデルの定義をする前に、lmfitのインポートは忘れないようにしてください。def関数でモデル関数を定義してから、model=lmf.Model(defで定義した関数名)のように書けば、modelという名前のモデル関数が定義できます。当然、modelという名前でなくは自分で別名を決めて問題ありません。
4. パラメータの初期値設定：モデル関数に対するパラメータを定義し、初期値を与えておきます。params=model.make_params(パラメータ1=数値1, パラメータ2=数値2,...)のように書くことで設定できます。paramsという名前は、自分で別名を決めて問題ありません。初期値は推定のものを入れれば良く、多くは直感や事前知識、視覚確認、簡単な計算により決定します。
5. フィッティングの実行：定義したモデルを用いて、用意したデータにフィットさせます。結果はresultに収納するため、result=model.fit(data, params, x)のようにして実行します。resultという名前も、自分で別名を決めて問題ありません。dataは用意したデータ（Y軸の値）、paramsは4で設定したパラメータ初期値、xは独立変数（X軸の値）です。独立変数xは基本、モデル関数の引数にすれば良いです。
6. フィッティング結果の表示：フィッティング結果を解析し、表示させます。print(result.fit_report())と書けば良いです。
7. データとフィット曲線のプロット：元のデータとフィット結果をプロットします。フィッティング曲線については、result.best_fitがフィッティング結果の配列になっており、プロットではY軸に対応します。

　たとえば、以下のように使います。ここではデータとして、 $y=3x+2$ の相関関係がある散布図を用意しています。この式から、モデル関数のパラメータが $a=3,b=2$ であることは明白ですが、いまは $a=1,b=1$ と適当に決めてフィッティングを行います。

#lmfitのインストール(省略)

#必要なモジュールのインポート
import numpy as np
import matplotlib.pyplot as plt
import lmfit as lmf

# データの準備
x = np.linspace(0, 10, 100)
y = 3.0 * x + 2.0 + np.random.normal(size=x.size)

# モデルの定義
def linear_model(x, a, b):
    return a * x + b

model = lmf.Model(linear_model)

# パラメータの初期値設定
params = model.make_params(a=1, b=1)

# フィッティングの実行
result = model.fit(y, params, x=x)

# フィッティング結果の表示
print(result.fit_report())

# データとフィット曲線のプロット
plt.scatter(x, y, label='Data')
plt.plot(x, result.best_fit, label='Best Fit', color='red')
plt.legend()
plt.show()

Ａ.２ リストの作り方

　データ解析ではリストの知識は不可欠です。以下では基本的な使い方を説明します。

Ａ.２.１ 空リスト

　要素数がゼロの空リストは、リストの名前を指定して、それに[]を代入すれば生成できます。リストの要素数はlen関数で取得できます。

list_empty = []

print(len(list_empty))
# 0

Ａ.２.２ 要素追加/削除

　要素の追加はappend関数、要素の削除はremove関数で可能です。以下では要素の追加・削除をするための関数を紹介します：

• 要素の追加

  1. append()
     　リストの末尾（最後）に要素を追加できます。文字列やリストもひとつの要素として追加され、結合はされません。

     l = [0, 1, 2]
     
     l.append(100)
     print(l)
     # [0, 1, 2, 100]
     
     l.append('abc')
     print(l)
     # [0, 1, 2, 100, 'abc']
     
     l.append([3, 4, 5])
     print(l)
     # [0, 1, 2, 100, 'abc', [3, 4, 5]]
     

  2. extend()
     　リストに別のリストやタプルを結合できます。すべての要素が元のリストの末尾に追加され、文字列やリストは各要素として一文字ずつ追加されます。
     　ちなみに、extend関数ではなく+演算子や+=演算子を使って連結することも可能です。+は新たなリストを返し、+=は既存のリストに追加します。

      l = [0, 1, 2]
     
     l.extend([10, 11, 12])
     print(l)
     # [0, 1, 2, 10, 11, 12]
     
     l.extend((100, 101, 102))
     print(l)
     # [0, 1, 2, 10, 11, 12, 100, 101, 102]
     
     l.extend('abc')
     print(l)
     # [0, 1, 2, 10, 11, 12, 100, 101, 102, 'a', 'b', 'c']
     
     l_new = l + [3, 4, 5]
     print(l_new)
     # [0, 1, 2, 10, 11, 12, 100, 101, 102, 'a', 'b', 'c', 3, 4, 5]
     
     l += [3, 4, 5]
     print(l)
     # [0, 1, 2, 10, 11, 12, 100, 101, 102, 'a', 'b', 'c', 3, 4, 5]
     

  3. insert()
     　リストの指定したインデックスの位置に要素を挿入できます。第一引数にインデックス、第二引数に挿入する要素を指定します。インデックスは、Pythonでは先頭から0, 1, 2, 3,...で、末尾から-1, -2, -3, -4,...と番号が振られますが、insert関数で挿入される位置は、-1だと末尾の一つ前となります。
     　また、append関数と同じく文字列やリストもひとつの要素として追加され、結合はされません。

     l = ['a', 'b', 'c']
     
     l.insert(1, 100)
     print(l)
     # ['a', 100, 'b', 'c']
     
     l.insert(0, 200)
     print(l)
     # [200, 'a', 100, 'b', 'c']
     
     l.insert(-1, 300)
     print(l)
     # [200, 'a', 100, 'b', 300, 'c']
     
     l.insert(-2, 400)
     print(l)
     # [200, 'a', 100, 'b', 400, 300, 'c']
     
     l.insert(0, [-1, -2, -3])
     print(l)
     # [[-1, -2, -3], 200, 'a', 100, 'b', 400, 300, 'c']
     

• 要素の削除
  1. clear()
     　リストのすべての要素が削除され、空のリストになります。

     l = [0, 1, 2]
     
     l.clear()
     print(l)
     # []
     

  2. pop()
     　リストの指定したインデックスの要素を削除し、さらに、その削除した要素の値を取得できます。insert関数でも述べましたが、インデックスは先頭から0, 1, 2, 3,...で、末尾から-1, -2, -3, -4,...と振られます。ちなみに、引数を省略して位置を指定しない場合は、末尾の要素が削除されます。

     l = [0, 10, 20, 30, 40, 50]
     
     popped_item = l.pop(0)
     print(popped_item)
     # 0
     
     print(l)
     # [10, 20, 30, 40, 50]
     
     popped_item = l.pop(3)
     print(popped_item)
     # 40
     
     print(l)
     # [10, 20, 30, 50]
     
     popped_item = l.pop(-2)
     print(popped_item)
     # 30
     
     print(l)
     # [10, 20, 50]
     

  3. remove()
     　リストの指定した値と同じ要素を検索し、最初の要素を削除できます。もし指定した値に一致する要素が複数含まれる場合には、最初の一つだけが削除されます。

     l = ['Alice', 'Bob', 'Charlie', 'Bob', 'Dave']
     
     l.remove('Alice')
     print(l)
     # ['Bob', 'Charlie', 'Bob', 'Dave']
     
     l.remove('Bob')
     print(l)
     # ['Charlie', 'Bob', 'Dave']
     

  4. del
     　リストの要素を範囲指定で削除できます。インデックスの指定はやはり先頭が0、末尾は-1です。範囲の指定はスライスで行い、複数の要素を一括で削除したり、全範囲を指定してすべての要素を削除したり、さらには範囲指定を[start:stop:step]として増分stepを指定すると、飛び飛びの複数の要素を一括で削除できます。

     l = [0, 10, 20, 30, 40, 50]
     
     del l[0]
     print(l)
     # [10, 20, 30, 40, 50]
     
     del l[3]
     print(l)
     # [10, 20, 30, 50]
     
     del l[-1]
     print(l)
     # [10, 20, 30]
     
     del l[-2]
     print(l)
     # [10, 30]
     
     ll = [0, 10, 20, 30, 40, 50]
     del ll[2:5]
     print(ll)
     # [0, 10, 50]
     
     ll = [0, 10, 20, 30, 40, 50]
     del ll[:3]
     print(ll)
     # [30, 40, 50]
     
     ll = [0, 10, 20, 30, 40, 50]
     del ll[-2:]
     print(ll)
     # [0, 10, 20, 30]
     
     lll = [0, 10, 20, 30, 40, 50]
     del lll[:]
     print(lll)
     # []
     
     llll = [0, 10, 20, 30, 40, 50]
     del llll[::2]
     print(llll)
     # [10, 30, 50]
     

• 配列の重ね合わせ
  1. numpy.hstack()
     　配列を横に重ねることができます。また、$n$ 次元配列の場合でも行数さえ同じであれば、重ねることが可能です。

     import numpy as np
     
     # np.arrange関数で公差1の等差数列を作成 
     arr1 = np.arange(3)
     print(arr1)
     # [0 1 2]
     
     arr2 = np.arange(3, 6)
     print(arr2)
     # [3 4 5]
     
     # np.hstack関数で配列を重ね合わせ
     arr3 = np.hstack((arr1, arr2))
     print(arr3)
     #[0 1 2 3 4 5]
     

  2. numpy.vstack()
     　配列を縦に重ねることができます。下の例で使っているnp.ones関数は、全要素を1で初期化した配列を生成するもので、np.ones((m, n))によって $m$ 行 $n$ 列の行列を指定できます。さらにnp.full関数は、全要素を任意の値で初期化した配列を生成するもので、np.ones((m, n),l)によって全要素がlとなった $m$ 行 $n$ 列の行列を指定できます。

     a1 = np.ones((2, 3), int)
     print(a1)
     # [[1 1 1]
     #  [1 1 1]]
     
     a2 = np.full((2, 3), 2)
     print(a2)
     # [[2 2 2]
     #  [2 2 2]]
     
     print(np.vstack([a1, a2]))
     # [[1 1 1]
     #  [1 1 1]
     #  [2 2 2]
     #  [2 2 2]]
     

Ａ.２.３ 要素番号

　enumerate関数を使うと、forループの中でリストやタプルなどのイテラブルオブジェクトの要素と同時に、インデックスを取得できます。インデックスは任意の値から開始可能です。

• for文
  　for 名前 in リストのようにすると、リストやタプルから要素をひとつずつ取り出すことができ、名前に格納されます。

  l = ['Alice', 'Bob', 'Charlie']
  
  for name in l:
      print(name)
  # Alice
  # Bob
  # Charlie
  

• enumerate関数
  　enumerate関数の引数にリストなどのイテラブルオブジェクトを指定すると、インデックス番号, 要素の順に取得できます。0以外の数値から開始したい場合は、enumerate関数の第二引数に任意の数値を指定すれば良いです。

  for i, name in enumerate(l):
      print(i, name)
  # 0 Alice
  # 1 Bob
  # 2 Charlie
  
  for i, name in enumerate(l, 1):
      print(i, name)
  # 1 Alice
  # 2 Bob
  # 3 Charlie
  

Ａ.２.４ 文字列の分割

　split関数を使うと、与えられた文字列に対して自分で設定したルールに基づく分割が可能となります。分割するには、リスト名.split('分割位置')のように書けば良いです。リスト名に入っている文字列を分割位置で分割して、結果を配列として格納することができます。さらに、リスト名.split('分割位置')[インデックス]とすれば、分割した文字列について、インデックスに対応する位置の文字列を引き出すこともできます。

file_list = ['/content/drive/MyDrive/IVdata_A5/IV_110mK_20230331.txt',
             '/content/drive/MyDrive/IVdata_A5/IV_115mK_20230331.txt',
             '/content/drive/MyDrive/IVdata_A5/IV_120mK_20230331.txt',
             '/content/drive/MyDrive/IVdata_A5/IV_125mK_20230331.txt',
             '/content/drive/MyDrive/IVdata_A5/IV_130mK_20230331.txt',
             '/content/drive/MyDrive/IVdata_A5/IV_135mK_20230331.txt',
             '/content/drive/MyDrive/IVdata_A5/IV_140mK_20230331.txt',
             '/content/drive/MyDrive/IVdata_A5/IV_145mK_20230331.txt',
             '/content/drive/MyDrive/IVdata_A5/IV_150mK_20230331.txt']
    
for tb in filename_list:
    tb_ = tb.split('/')[-1]
    print(tb_)
    '''
    IV_110mK_20230331.txt
    IV_115mK_20230331.txt
    IV_120mK_20230331.txt
    IV_125mK_20230331.txt
    IV_130mK_20230331.txt
    IV_135mK_20230331.txt
    IV_140mK_20230331.txt
    IV_145mK_20230331.txt
    IV_150mK_20230331.txt
    '''
  for i, name in enumerate(tb_.split('_')):
      print(i, name)
      # 下の出力例はIV_110mK_20230331.txt のforループのみ
      '''
      1 IV
      2 110mK
      3 20230331.txt
      '''

Ａ.２.５ 文字列の取得

　find関数を使うと、文字列中にある特定の文字の位置を取得できます。find関数の第一引数に指定した文字が文字列に含まれている場合は、その文字に対する最初の位置が返されます。また、指定した文字が含まれていない場合には、-1が返されます。
　文字列にスペースが含まれている場合、スペースも一文字としてカウントされることに注意してください。

TES = 'Transition Edge sensor'

print(TES.find('Sensor'))
# 17

print(TES.find('X-ray'))
# -1

Ａ.２.６ 熱浴温度リスト

　さて、今回定義したmakeTbath関数で行っているプログラムを下の図にまとめました。図中では、具体的に一つ目のファイルのスキャンで行っていることをまとめてあります。

def makeTbath(filename_list):
  Tbath = []
  num = 0
  for tb in filename_list:
    tb_ = tb.split('/')[-1]
    for e, name in enumerate(tb_.split('_')):
      if name.find('mK')!=-1:
        num = e

    tb_key = tb_.split('_')[num]
    tbath = float(tb_key[:-2])
    Tbath.append(tbath)

  Tbath = np.asarray(Tbath)
  print(Tbath)
  return Tbath