Chapter 16

플라즈마 모델링: 유체·PIC·하이브리드

8장 전역 모델은 챔버 전체를 숫자 몇 개(Tₑ, n₀)로 줄였다. 밀도 분포, 시스 구조, 웨이퍼 위 균일도, 이온 에너지 분포를 알려면 공간과 속도 공간을 나눠 풀어야 한다. 전자를 연속체(유체)로 보면 빠르지만 분포 함수의 꼬리를 잃고, 입자로 추적하면(PIC/MCC) 정확하지만 데바이 길이와 플라즈마 주기를 다 풀어야 해서 비싸다. 이 장은 표류-확산 유체 풀이와 1D PIC를 직접 돌려 보며 각 모델이 무엇을 예측하고 언제 깨지는지, 수치 조건을 어기면 어떤 가짜 물리가 생기는지, 그리고 장비 설계에 쓰는 하이브리드 모델이 어떻게 조립되는지 본다.

전역·유체·운동론 모델의 가정, 출력, 계산 비용을 비교하고 문제에 맞는 모델을 고를 수 있다.
표류-확산 근사와 포아송 방정식을 결합한 1D RF 방전 모델을 세우고, 샤페터-거멜 이산화와 반암시적 포아송이 필요한 이유를 설명한다.
국소장 근사(LFA)와 국소 평균 에너지 근사(LMEA)의 차이를 음극 시스·글로에서 이온화 원천으로 설명한다.
PIC/MCC 한 사이클(전하 할당 → 장 풀이 → 보간 → 입자 밀기 → 충돌)을 따라가고, 안정 조건 $\Delta x \lesssim \lambda_D$ , $\omega_{pe}\Delta t \lesssim 0.2$ 을 어겼을 때의 수치 가열·발산을 알아본다.
ICP 하이브리드 모델(전자기 → 전자 에너지 → 유체 수송)의 반복 구조와, 모델을 진단으로 검증·교정하는 절차와 함정을 설명한다.

모델의 위계

플라즈마의 모든 정보는 종마다 위치·속도의 분포 함수 $f_s(\mathbf r, \mathbf v, t)$ 에 들어 있고, 이것이 볼츠만 방정식을 따른다. 모델은 이 방정식을 얼마나 줄이느냐로 나뉜다.

\frac{\partial f_s}{\partial t} + \mathbf v\cdot\nabla_{\mathbf r}f_s + \frac{q_s}{m_s}(\mathbf E + \mathbf v\times\mathbf B)\cdot\nabla_{\mathbf v}f_s = \left(\frac{\partial f_s}{\partial t}\right)_{c}

6차원(+시간) 방정식. 운동론(PIC/MCC, 직접 볼츠만 풀이)은 이것을 거의 그대로 푼다. 유체는 속도에 대해 적분한 모멘트(밀도, 플럭스, 에너지)만 풀고 분포 모양은 가정한다. 전역은 공간까지 평균해 상미분 방정식 몇 개로 줄인다.

그림 16-1. 모델 위계. 아래로 갈수록 빠르고 위로 갈수록 정확하다. 실무에서는 전역 모델로 조건을 훑고, 유체·하이브리드로 장비 형상을 설계하고, PIC/MCC로 핵심 물리(저압 시스, IEDF)를 확인한다.

모델 선택의 첫 잣대는 크누센 수 $\lambda/L$ (평균 자유 행로 / 관심 길이)다. 전자 에너지 이완 길이가 장치보다 짧으면(수 Torr PECVD) 전자가 국소 전기장과 평형을 이뤄 유체가 잘 맞는다. 수 mTorr ICP·CCP에서는 전자 에너지 이완 길이가 챔버만 해 EEDF가 비국소적이고(4장), 시스를 건너는 이온도 충돌이 거의 없어 IEDF를 유체로 얻을 수 없다.

표류-확산 유체 모델과 포아송

운동량 방정식에서 관성을 버리고(충돌이 잦아 가속 시간 ≪ 관심 시간) 압력 기울기와 전기력이 충돌 마찰과 균형을 이룬다고 보면 플럭스가 표류-확산 근사(drift-diffusion approximation)로 닫힌다. 전자·이온 연속 방정식과 포아송 방정식이 서로를 물고 돈다.

\frac{\partial n_s}{\partial t} + \nabla\cdot\boldsymbol\Gamma_s = S,\qquad \boldsymbol\Gamma_e = -\mu_en_e\mathbf E - D_e\nabla n_e,\qquad \boldsymbol\Gamma_i = \mu_in_i\mathbf E - D_i\nabla n_i,\qquad \nabla\cdot(\varepsilon_0\nabla\phi) = -e(n_i - n_e)

\mathbf E = -\nabla\phi

. 국소장 근사에서 이온화 원천은

S = \alpha(E/p)\,|\boldsymbol\Gamma_e|

\alpha = Ap\,e^{-Bp/|E|}

(Ar:

A \approx 12

cm⁻¹Torr⁻¹,

B \approx 180

V cm⁻¹Torr⁻¹). 이동도

\mu_ep \approx 30

m²Torr/(V s),

\mu_ip \approx 0.14

m²Torr/(V s), 확산 계수는 아인슈타인 관계

D = \mu T

(6장).

그림 16-2. 표류-확산 모델의 한 시간 단계. 장이 플럭스를, 플럭스가 밀도를, 밀도가 다시 장을 정한다. 전자는 가볍고 빨라 이 고리가 강하게 결합되어 있다. 순서대로 따로 풀면(분리) 시간 간격이 유전 이완 시간에 묶인다.

수치로 풀 때 두 가지 함정이 있다. 첫째, 유전 이완 시간 $\tau_d = \varepsilon_0/(e\mu_en_e)$ 은 1 Torr에서 10¹⁵ m⁻³일 때 약 2 ns, 10¹⁷ m⁻³일 때 약 20 ps로 밀도에 반비례한다. 포아송을 앞 단계 밀도로 풀면(명시적) 시간 간격이 $\tau_d$ 보다 작아야 하고, 아니면 전하가 출렁이며 발산한다. 다음 단계의 표류 전류를 포아송에 미리 넣는 반암시적 방법이 이 제한을 없앤다. 둘째, 시스 안에서는 표류가 확산보다 수십 배 커서(셀 페클레 수 $|a| = |v|h/D \gg 2$ ) 중앙 차분이 음의 밀도와 진동을 만든다. 반도체 소자 시뮬레이션에서 온 샤페터-거멜(Scharfetter–Gummel) 이산화가 표준 해법이다(TCADBook의 소자 방정식과 같은 꼴).

\Gamma_{i+1/2} = \frac{D}{h}\Big[B(-a)\,n_i - B(a)\,n_{i+1}\Big],\qquad B(x) = \frac{x}{e^x - 1},\qquad a = \frac{v\,h}{D};\qquad \nabla\cdot\!\big[(\varepsilon_0 + e\Delta t(\mu_en_e + \mu_in_i))\nabla\phi^{k+1}\big] = -\rho^{k}

셀 안에서 플럭스가 일정하다고 보고 표류-확산 방정식을 정확히 푼 지수 꼴의 해를 이은 것이다.

a \to 0

이면 중앙 차분 확산,

|a| \to \infty

이면 바람 거슬러(upwind) 표류로 자연스럽게 넘어간다. 암시적 연속 방정식과 함께 쓰면 계수 행렬이 M-행렬이 되어 밀도가 음이 되지 않는다.

그림 16-3. 샤페터-거멜은 격자점 사이 밀도를 직선이 아닌 표류-확산 방정식의 정확한 지수 해로 잇는다. 그래서 큰 페클레 수에서도 진동 없이 표류 플럭스를 옮긴다.

SIMULATOR

이산화 비교: 중앙 차분 · 바람 거슬러 · 샤페터-거멜

전역 페클레 수 Pe = vL/D 셀 수 N

셀 페클레 수 a—

최대 오차: 중앙—

최대 오차: 바람 거슬러—

최대 오차: SG—

정상 1D 표류-확산

d(vn - D\,dn/dx)/dx = 0

n(0) = 1

n(1) = 0

의 정확한 해

n = (e^{Pe} - e^{Pe\,x})/(e^{Pe} - 1)

(회색)와 세 이산화의 격자 해(점). 해볼 것: ① N = 8, Pe = 20(a = 2.5)에서 중앙 차분은 음의 밀도로 출렁인다(빨강). ② 바람 거슬러는 늘 양수지만 수치 확산으로 경계층이 퍼진다. ③ SG는 셀 수와 상관없이 격자점에서 정확하다. 실제 시스는 a가 10–100이다.

1D RF 방전을 직접 풀기

위 식들을 평행판 사이에서 풀어 보자. 왼쪽 전극에 13.56 MHz 전압 $V_{rf}\sin\omega t$ , 오른쪽은 접지, 두 전극에서 밀도 0(흡수). 시간 간격은 RF 주기의 1/100(약 0.74 ns)로, 이 모델이 도달하는 밀도(수×10¹⁵ m⁻³)의 유전 이완 시간(약 0.5 ns)보다 길고 방전이 더 조밀해지면 수십 배 길어지지만 반암시적 포아송과 암시적 SG 덕분에 안정하다. 처음에 작은 씨앗 밀도에서 시작하면 전자 사태로 밀도가 자라고, 전극 앞에서 전자가 빠져 양전하층(시스)이 생기며, 시스가 벌크 전기장을 가려 이온화가 줄어 정상 상태에 이른다. 이 모델은 수백 mTorr–수 Torr(PL-300 PE 모듈 PECVD 조건) 쪽에서 믿을 만하고, 수십 mTorr 아래로 내려가면 국소장 근사가 깨진다.

SIMULATOR

1D 표류-확산 + 포아송: RF 방전이 시스를 만드는 과정

전극 간격

압력 RF 전압 진폭 프레임당 RF 주기

흐른 RF 주기—

최대 nᵢ(주기 평균)—

시스 두께(왼쪽)—

전극 이온 플럭스—

상태—

Ar, Tₑ = 2 eV(확산 계수용 고정), 101 격자점, 주기당 100 시간 단계, 국소장 이온화

\alpha(E/p)|\Gamma_e|

. 위: 주기 평균 nₑ(파랑 계열)·nᵢ(이온 색), 얇은 선은 현재 순간의 nₑ. 아래: 현재 위상의 전위 φ(x)(진한 선)와 한 주기 동안의 포락선(옅은 띠). 해볼 것: ① 처음 수십 주기 동안 밀도가 지수 함수로 자라다 시스가 생기며 포화한다. ② 압력을 3 Torr로 올리면 시스가 얇아지고 밀도 봉우리가 뾰족해진다. ③ 전압을 250 → 100 V로 낮추면 밀도가 약 1/3로 준다. 이 모델에는 2차 전자 방출과 전자 에너지 수송이 없어 실제 유지 전압 문턱은 잘 맞추지 못한다. ④ 10 mm 간격(PECVD)에서는 시스가 간격의 큰 몫을 차지한다. 시스 두께는 주기 평균 nₑ가 nᵢ의 절반이 되는 위치로 잰다.

국소장 근사와 국소 평균 에너지 근사

위 모델은 이온화율을 그 자리의 전기장 $E/p$ 하나로 정했다. 이것이 국소장 근사(local field approximation, LFA)다. 전자가 장과 평형을 이루는 길이(에너지 이완 길이 $\lambda_\varepsilon$ )가 장의 변화 길이보다 짧을 때만 맞다. 음극 시스처럼 장이 수 mm 안에서 0으로 떨어지는 곳에서는 시스에서 에너지를 얻은 전자가 장이 거의 없는 글로 영역까지 날아가 이온화한다. LFA는 이 이온화를 놓치고(글로에서 원천 0), 시스 깊숙한 곳의 원천을 과대평가한다.

국소 평균 에너지 근사(local mean energy approximation, LMEA)는 전자 평균 에너지 $\bar\varepsilon$ 에 대한 수송 방정식을 하나 더 풀고, 반응 속도와 수송 계수를 $E/p$ 가 아닌 $\bar\varepsilon$ 의 함수로 표에서 찾는다(표는 볼츠만 풀이기로 미리 만든다).

\frac{\partial(n_e\bar\varepsilon)}{\partial t} + \nabla\cdot\boldsymbol\Gamma_\varepsilon = -e\boldsymbol\Gamma_e\cdot\mathbf E - n_eN\,\mathcal L(\bar\varepsilon)\quad\Rightarrow\quad \frac{d\bar\varepsilon}{dx} = E(x) - E_{eq}(\bar\varepsilon)\ \ (\text{1D 정상 빔 근사})

\mathcal L

은 탄성·비탄성 충돌로 잃는 에너지 손실 함수,

E_{eq}(\bar\varepsilon)

는 그 평균 에너지를 국소 평형으로 유지하는 장(여기서는

E_{eq} = c\,p\,\bar\varepsilon^{3/2}

c \approx 316

V m⁻¹Torr⁻¹eV^−3/2로 근사). 장이 갑자기 떨어지면

\bar\varepsilon

는 이완 길이

\lambda_\varepsilon \sim \bar\varepsilon/E_{eq}

에 걸쳐 천천히 식는다.

그림 16-4. 음극 시스와 글로. LFA는 이온화 원천을 장이 큰 시스 안에 가두고 글로에서 0으로 만든다. 에너지 수송을 푸는 LMEA(또는 운동론 모델)에서는 시스에서 가속된 전자가 글로로 들어가 이온화가 번진다. 실제 음 글로가 빛나는 이유다.

SIMULATOR

LFA vs LMEA: 음극 시스 뒤 이온화 원천

압력 음극 강하 V_c 시스 두께 d_s

이완 길이 λ_ε(10 eV)—

원천 최대 위치 LFA · LMEA—

글로 쪽 이온화 몫 LFA · LMEA—

증배 ln M LFA · LMEA—

선형 시스 장

E = (2V_c/d_s)(1 - x/d_s)

, 글로 장

E_b = 200p

V/m. 전자는 음극에서 1 eV로 출발해

d\bar\varepsilon/dx = E - cp\bar\varepsilon^{3/2}

를 따르고, 이온화 계수는 LFA에서

\alpha(E)

, LMEA에서

\alpha(E_{eq}(\bar\varepsilon))

다. 원천

S = \alpha\Gamma

d\Gamma/dx = \alpha\Gamma

. 위: 정규화한 원천 S(x)(LFA 점선, LMEA 실선)와 시스 경계. 아래: 평균 에너지 ε̄(x)와 국소 평형 에너지 ε_loc(E(x)). 해볼 것: ① 0.2 Torr로 내리면 λ_ε가 시스보다 길어져 LMEA 원천이 글로 깊숙이 들어간다. ② 5 Torr에서는 λ_ε ≪ d_s라 두 근사의 원천 봉우리 위치가 거의 같아진다. ③ LFA는 음극 바로 앞(x ≈ 0)이 장이 가장 커서 α도 가장 크다고 보지만, 실제로는 막 출발한 전자의 에너지가 낮아 거기서 이온화하지 못한다(아래 그림에서 ε̄가 ε_loc를 늦게 따라감).

PIC/MCC: 입자로 풀기

PIC(particle-in-cell)은 실제 전자·이온 수천억 개를 대신하는 거대 입자(macro-particle, 한 개가 10⁶–10⁹개를 대표) 수만–수억 개를 추적한다. 입자는 연속 공간에서 움직이지만 전기장은 격자에서 푼다. 충돌은 MCC(Monte Carlo collision)로 확률적으로 처리한다. 한 시간 단계는 네 단계(+충돌)의 고리다.

\rho_j = \sum_p q_p\,W(x_p - x_j),\quad \frac{\phi_{j+1} - 2\phi_j + \phi_{j-1}}{\Delta x^2} = -\frac{\rho_j}{\varepsilon_0},\quad v_p^{n+1/2} = v_p^{n-1/2} + \frac{q}{m}E(x_p^n)\Delta t,\quad x_p^{n+1} = x_p^n + v_p^{n+1/2}\Delta t

W

는 할당 함수(구름 입자법 CIC: 이웃 두 격자점에 거리에 따라 선형으로 나눔, 가까운 점이 더 많이 받음). 장을 입자 위치로 보간할 때도 같은

W

를 써야 자기 힘이 사라진다. 속도와 위치를 반 단계 엇갈려 두는 립프로그(leapfrog)는 2차 정확도이고 에너지를 장기간 잘 보존한다. 충돌 확률은

P = 1 - e^{-\nu\Delta t}

이고, 에너지마다 다른

\nu(\varepsilon)

는 최대 주파수

\nu_{max}

로 뽑은 뒤

\nu(\varepsilon)/\nu_{max}

확률로 받아들이는 널 충돌법(null collision)으로 효율을 높인다.

그림 16-5. PIC/MCC 사이클. 입자 → 격자(할당), 격자에서 장 풀이, 격자 → 입자(보간), 입자 밀기, 그리고 충돌. 계산 시간의 대부분은 ①③④처럼 입자 수에 비례하는 단계가 차지한다.

SIMULATOR

PIC 한 사이클 따라가기(입자 12개, 격자 8칸)

자동 재생 충돌 주파수 ν/ω_pe

단계—

시간 tω_pe—

전체 에너지—

주기 경계, 정규화 단위(ω_pe = 1, 이온은 균일한 양전하 바탕), Δt = 0.3. 위: 입자(점, 화살표 = 속도)와 격자점, 막대 = 할당된 전하 밀도 ρ_j. 아래: 전위 φ_j(선)와 전기장 E_j(점선). 한 번 누를 때마다 할당 → 장 풀이 → 보간(입자 위 화살표 = 힘) → 밀기 → 충돌(주황 테두리 = 이번에 충돌한 입자, 속도 방향이 무작위로 바뀜) 순서로 진행한다.

안정 조건: 무엇을 풀어야 하나

PIC의 수치 조건은 물리 척도에서 나온다(2장). 격자가 데바이 길이를 못 풀면( $\Delta x \gg \lambda_D$ ) 격자보다 짧은 파장의 요동이 긴 파장으로 둔갑(에일리어싱)해 전자를 데우는 유한 격자 불안정이 생기고, $\lambda_D$ 가 $\Delta x/\pi$ 정도로 커질 때까지 가짜 가열이 이어진다. 시간 간격은 립프로그가 플라즈마 진동을 풀어야 하므로 $\omega_{pe}\Delta t \lt 2$ 가 안정 한계, 정확도를 위해 $\omega_{pe}\Delta t \lesssim 0.2$ 다. 입자가 한 단계에 한 칸 이상 건너뛰면 안 되므로( $v\Delta t \lt \Delta x$ ) 고에너지 전자도 고려한다. 셀당 입자 수가 적으면 통계 잡음이 $1/\sqrt{N_{ppc}}$ 로 커져 충돌처럼 작용하는 수치 충돌·가열이 생긴다.

\Delta x \lesssim \lambda_D = \sqrt{\frac{\varepsilon_0T_e}{en_e}},\qquad \omega_{pe}\Delta t \lesssim 0.2\ \ (\text{안정 한계 } 2),\qquad \sin\frac{\omega\Delta t}{2} = \frac{\omega_{pe}\Delta t}{2},\qquad \omega^2 = \omega_{pe}^2 + 3k^2v_{th}^2

세 번째 식은 립프로그가 차가운 플라즈마 진동을 풀 때의 수치 분산 관계다.

\omega_{pe}\Delta t

가 2를 넘으면 실수해가 없어 진폭이 매 단계 커진다. 마지막은 따뜻한 플라즈마의 보옴-그로스 분산 관계다.

그림 16-6. PIC의 두 조건. 격자는 데바이 길이를, 시간 간격은 플라즈마 진동을 풀어야 한다. PL-300 IC 기준 플라즈마에서

\lambda_D \approx 19

µm,

1/\omega_{pe} \approx 28

ps이므로 명시적 PIC로 챔버 전체(0.25 m)를 ms 동안 푸는 것은 사실상 불가능하다.

SIMULATOR

1D 정전 PIC: 플라즈마 진동, 이류 불안정, 수치 가열

실험

격자 Δx/λ_D 시간 간격 ω_pe Δt 셀당 입자 수 충돌(MCC) ν/ω_pe

시간 tω_pe—

전체 에너지 변화—

측정 ω/ω_pe (이론)—

Tₑ/Tₑ₀ · 성장률 γ—

64칸 주기 경계, 전자 거대 입자 64 × (셀당 입자 수), 이온은 균일 바탕, CIC 할당·보간, 립프로그. 정규화: ω_pe = 1, λ_D = v_th/ω_pe = 1. 위: 위상 공간(x, v). 아래: 운동 에너지(전자 색)·장 에너지(강조색, 로그)·전체(회색), 처음 전체 에너지로 정규화. 해볼 것: ① 플라즈마 진동은 장 에너지가 2ω로 출렁이며, 측정 ω는 보옴-그로스 값(립프로그 보정 포함)과 맞는다. ω_peΔt를 1.5로 키우면 주파수가 눈에 띄게 어긋나고, 2를 넘기면 발산한다. ② 이류 불안정은 장 에너지가 로그 그래프에서 직선(γ ≈ 0.35 ω_pe)으로 자라다 위상 공간에 소용돌이가 생기며 포화한다. ③ 열 플라즈마에서 Δx/λ_D를 3 이상으로 두면 Tₑ가 서서히 올라간다(유한 격자 가열). 0.5로 두면 거의 일정하다. ④ 충돌을 넣으면 진동이 감쇠한다.

하이브리드와 2D 장비 모델

장비 설계용 2D·3D 모델은 대부분 하이브리드다. 서로 시간 척도가 다른 물리를 모듈로 나누고, 각 모듈을 자기 척도로 풀어 결과를 주고받으며 수렴할 때까지 반복한다. ICP 모델의 전형적인 구성은 다음과 같다.

전자기 모듈: 코일 전류가 만드는 방위각 장을 주파수 영역 헬름홀츠 식 $\nabla^2E_\theta - E_\theta/r^2 + k_0^2\varepsilon_p E_\theta = i\omega\mu_0J_{coil}$ 로 푼다. 플라즈마 유전율 $\varepsilon_p = 1 - \omega_{pe}^2/(\omega(\omega - i\nu_m))$ 이 밀도에 달려 있어 표피 깊이(10장)가 결과로 나온다.
전자 에너지 모듈: 흡수 전력 분포 $\tfrac12\mathrm{Re}(\sigma)|E|^2$ 로 전자 온도(유체 에너지 방정식) 또는 EEDF(전자 몬테카를로)를 구하고, 이온화·해리 속도를 낸다.
유체 수송 모듈: 이온·중성·래디컬의 연속·운동량 방정식, 양극성 확산, 기체 흐름과 펌핑. 전자는 준중성으로 묶는다. 이온은 시스 입구까지 유체로, 시스 안과 웨이퍼 IEDF는 별도 이온 몬테카를로로 뽑는다.
회로·화학 모듈: 정합기·블로킹 커패시터·바이어스 회로, 수십–수백 개 반응 화학과 표면 반응(벽 재결합, 식각 생성물).

그림 16-7. ICP 하이브리드 모델의 모듈 고리. 전자기 모듈의 흡수 전력이 이온화 원천을, 원천이 밀도를, 밀도가 다시 유전율과 표피 깊이를 정한다. 각 모듈은 자기 시간 척도(ns~ms)로 따로 풀고 결과만 주고받는다.

아래 시뮬레이터는 이 고리를 가장 단순하게 줄인 것이다. PL-300 IC 모듈의 (r, z) 단면에서 ① 코일 아래 전력 흡수를 표피 깊이 $\delta = c/\omega_{pe}$ (무충돌 근사)로 감쇠시키고 ② 이온화 원천 $S = P/(e\mathcal E_T)$ 로 양극성 확산 $\nabla\cdot(D_a\nabla n) + S = 0$ 을 풀고(벽에서 보옴 플럭스 경계) ③ 새 밀도로 표피 깊이를 다시 구하는 반복을 수렴할 때까지 돈다.

그림 16-8. 축대칭 2D 계산 영역. 반쪽만 풀고 거울상으로 보인다. 위는 세라믹 창(코일), 아래는 웨이퍼, 옆은 벽. 상용 코드는 코일·창·척·펌핑 포트까지 수만–수십만 개의 비정렬 격자로 나눈다.

SIMULATOR

2D ICP 하이브리드 반복: 표피 흡수 ↔ 양극성 확산

흡수 전력 압력 바깥 코일 전류 비중

반복 · 변화량—

최대 nₑ—

표피 깊이 δ(중간 반지름)—

웨이퍼 이온 플럭스 균일도—

D_a

는 Ar⁺ 이동도와 Tₑ = 2.5 eV로

D_a \approx 38\times(10\,\text{mTorr}/p)

m²/s,

\mathcal E_T

≈ 108 eV(Tₑ = 2.5 eV, 전역 모델 값), 벽 경계 플럭스

n/(\Delta/2D_a + 1/u_B)

. 코일은 반지름 8 cm(안쪽)·18 cm(바깥) 두 고리로, 전력 흡수는

e^{-2(H-z)/\delta(r)}

로 감쇠, δ(r)는 창 아래 밀도로 구한다(반복마다 반만 갱신). 위: n(r, z) 단면(좌우 거울상, 흰 선 = 흡수 전력 밀도 50 % 등고선). 아래: 웨이퍼 이온 플럭스 Γ(r)(웨이퍼 영역 강조)와 반복에 따른 최대 밀도 수렴. 해볼 것: ① 전력을 올리면 밀도가 커지고 δ가 얇아져 흡수가 창에 붙는다. ② 바깥 코일 비중으로 균일도를 맞춰 보라. 최선은 약 0.45–0.5에서 ±4 % 남짓으로, 두 고리만으로는 PL-300 목표 ±3 %에 닿지 않는다. 실제 장비가 코일을 3–4 구역으로 나누고 구역별 전류를 따로 조절하는 이유다. ③ 압력을 올리면 확산이 느려져 분포가 원천 모양을 더 닮는다. 이 단순 모델의 중심 밀도(1.5 kW, 10 mTorr에서 약 2.4×10¹¹ cm⁻³)는 PL-300 전역 모델 값(4×10¹¹)보다 약 40 % 낮다. 저압에서 확산 근사와 보옴 경계가 벽 손실을 과대평가하기 때문이며, 이런 차이를 메우는 것이 다음 절의 교정이다.

검증과 교정: 모델을 진단과 맞대기

모델의 입력에는 불확실한 값이 많다(단면적, 벽 재결합 계수, 2차 전자 계수, 전력 전달 효율). 교정(calibration)은 그런 손잡이를 측정에 맞게 조정하는 것이고, 검증(validation)은 교정에 쓰지 않은 다른 조건·다른 진단으로 예측을 시험하는 것이다. 둘을 섞으면 "맞춘 것을 맞혔다"는 자기만족이 된다. 진단 자체의 오차(15장: 탐침 n은 ±30 % 수준)도 비교에 넣어야 한다. 탐침의 계통 오차가 모델 매개변수에 흡수되면 다른 조건에서 예측이 빗나간다.

그림 16-9. 교정과 검증. 교정에 쓴 자료로 모델을 평가하지 않는다. 가능하면 원리가 다른 진단(탐침 vs 간섭계, OES vs 질량 분석)으로 교차 확인한다.

SIMULATOR

모델 교정: 탐침 지름 분포에 확산 모델 맞추기

모델: 원천 반지름 r_s 모델: 확산 계수 배율 탐침 계통 오차 +25 %(해석 이론)

χ²/자유도—

r_s · D 배율—

검증 조건 중심 n 오차—

모델: 1D 지름 방향 정상 확산

(1/r)(rD\,n')' + S(r) = 0

, 가우스 고리 원천(폭 4 cm, 전력 1.5 kW 고정), 벽에서 보옴 플럭스. "측정": 숨겨진 참 매개변수(r_s = 0.13 m, D 배율 1.4)로 만든 분포에 3 % 잡음을 더한 탐침 13점, 오차 막대 ±5 %. 검증은 압력을 3배(D를 1/3)로 바꾼 조건의 중심 밀도를 교정된 모델로 예측해 참값과 비교한다. 해볼 것: ① 손으로 χ²/자유도를 1 근처까지 낮춰 보라. ② 자동 교정은 참값 근처를 찾는다. ③ 계통 오차를 켜고 교정하면 χ²는 여전히 좋지만 D 배율이 틀린 값(약 1.1)으로 가고, 검증 조건에서도 25 % 틀린다. 측정과 잘 맞는 모델이 옳은 모델이라는 보장은 없다.

계산 비용과 모델 고르기

명시적 PIC의 비용은 격자 수 × 셀당 입자 수 × 시간 단계 수에 비례하고, 격자와 시간 간격이 $\lambda_D$ 와 $1/\omega_{pe}$ 에 묶여 있어 밀도가 오르면 빠르게 커진다(1D에서 $\propto n^{1/2}\cdot n^{1/2}$ , 3D에서 $\propto n^{3/2}\cdot n^{1/2} = n^2$ ). 그래서 실제 PIC는 밀도를 낮추거나 장치를 줄이고(축척), 질량비를 바꾸고, 암시적·에너지 보존 PIC나 GPU를 쓴다. 유체는 시스를 풀 만큼만 격자를 촘촘히 하면 되고 반암시적으로 큰 시간 간격을 쓸 수 있어 수 자릿수 싸다.

\text{비용} \approx N_{cell}\,N_{ppc}\,N_{step}\,c_p,\qquad N_{cell} = \prod_d\frac{L_d}{\Delta x},\qquad N_{step} = \frac{t_{sim}}{\Delta t} = \frac{t_{sim}\,\omega_{pe}}{0.2}

c_p

는 입자 하나를 한 단계 처리하는 연산 수(할당·보간·밀기·충돌 약 50–100 flop). 유체는

N_{cell}N_{step}c_f

c_f

는 셀 하나에 방정식 여럿을 암시적으로 푸는 비용(수천 flop).

SIMULATOR

계산 비용 추정기: 격자 × 입자 × 시간 단계

차원

전자 밀도 모사 시간 PIC 셀당 입자 수 코어 수 밀도 축척(PIC)

λ_D · 1/ω_pe—

PIC 격자 · 입자—

PIC 시간 단계—

PIC 벽시계 시간—

영역: PL-300 IC 모듈(1D: 창-웨이퍼 0.10 m, 2D: 반지름 0.25 × 높이 0.10 m, 3D: 0.5 × 0.5 × 0.10 m 상자), Tₑ = 2.5 eV. PIC: Δx = λ_D, ω_peΔt = 0.2, 전자·이온 두 종, 입자-단계당 80 flop, 장 풀이 셀-단계당 30 flop. 유체: Δx = max(λ_D, 0.5 mm)(시스 근처를 풀 정도), Δt = RF 주기/200, 셀-단계당 2000 flop. 하이브리드: 유체 이온(Δt = 1 ns, 셀-단계당 1000 flop) + 전자 몬테카를로 10⁵개 × 10⁴ 단계 × 반복 50회. 전역: 약 10⁸ flop. 코어당 유효 1 GFLOP/s, 병렬 효율 70 %. "밀도 축척"은 PIC에서 밀도를 10^s 배 낮춰 푸는 흔한 꼼수(λ_D ↑, ω_pe ↓). 막대는 벽시계 시간(로그). 해볼 것: ① 기준 플라즈마의 2D PIC는 100 µs만 풀어도 64코어로 십 년이 넘게 걸린다. 밀도 축척 2–3을 넣어야 일 단위가 된다(대신 시스·EEDF 물리가 바뀐다). ② 1D PIC는 노트북으로도 충분하다. 저압 CCP의 IEDF·EEDF 연구가 1D PIC로 많이 이뤄진 이유다.

질문	권장 모델	이유
압력·전력·기체 비율을 넓게 훑어 Tₑ, n, 래디컬 경향 보기	전역 모델	초 단위, 수백 반응 화학도 가능
수 Torr PECVD 평행판의 시스·밀도 분포	1D/2D 유체(LFA 또는 LMEA)	λ_ε ≪ 간격, 국소 평형이 맞음
ICP 코일·창 설계와 웨이퍼 균일도	2D 하이브리드(EM + 유체 + 전자 에너지)	형상·전자기·화학이 함께 필요
저압 CCP의 IEDF, 확률 가열, 이중 주파수 결합	1D PIC/MCC	무충돌 시스·비국소 EEDF는 운동론만 정확
VHF 정상파·표피 효과 균일도	2D/3D 전자기 + 유체	파장이 전극만 해 전자기 풀이가 핵심
고종횡비 구멍 안 이온·전자 그늘, 대전	시스 IEDF(PIC/MC) + 형상 내 입자 추적	구멍 척도 µm, 반응기 모델과 연결(TCADBook 형상 모사)
펄스 잔광 µs 동역학	시간 의존 전역 + 1D PIC	Tₑ 급랭은 운동론, 밀도 감쇠는 전역으로 충분

조건 변화	모델에 미치는 결과
압력 ↓ (1 Torr → 10 mTorr)	λ_ε·λᵢ가 장치만 해짐 → LFA 붕괴, 유체 IEDF 무의미, PIC 필요성 ↑. 대신 밀도가 낮아 PIC 비용은 덜 큼
전력 ↑ (nₑ 10배)	λ_D가 1/√10, ω_pe가 √10배 → 2D PIC 비용 약 30–100배, 유체는 시스 격자만 조금 촘촘히
주파수 ↑ (13.56 → 60 MHz)	유체·PIC 모두 RF 주기 해상 비용 ↑, 전자기(정상파) 모듈 필요
전기 음성 기체 추가	음이온 종과 반응 수십 개 추가, 이온-이온 플라즈마는 준중성 가정이 바뀜 → 유체·전역 모두 화학이 비용 지배
격자 Δx를 λ_D의 3배로 키움(PIC)	유한 격자 가열로 Tₑ가 가짜로 오름 → 계산이 빨라진 대가로 물리가 틀림

IC 모듈 균일도 개발의 모델 사슬코일 배치 후보 수십 개는 2D 하이브리드 모델로 하루 안에 거르고(시뮬레이터처럼 바깥 코일 비중을 훑어 웨이퍼 이온 플럭스 ±3 % 안의 창을 찾는다), 상위 3개를 3D로 확인한다(코일 끝단·급전부가 만드는 방위각 비대칭, 2D로는 볼 수 없다). 바이어스 IEDF는 시스만 떼어 1D PIC/MCC로 따로 푼다. 마지막에 웨이퍼형 RFEA와 탐침 스캔(15장)으로 검증하는데, 탐침 절대값(±30 %) 대신 분포 모양과 전력에 따른 기울기를 비교한다.

PIC로 PL-300을 바로 풀 수 없는 이유기준 플라즈마

\lambda_D \approx 19

µm,

\omega_{pe} \approx 3.6\times10^{10}

s⁻¹. 2D (0.25 × 0.10 m) 명시적 PIC는 격자 약 7×10⁷개, 셀당 100입자면 입자 약 10¹⁰개, 1 ms를 풀려면 약 1.8×10⁸ 단계다. 연산 약 10²⁰ flop로, 1000 코어에서도 10년 가까이 걸린다. 실제 연구는 1D·축척·암시적 PIC로 한다.

핵심 정리

모델 위계: 전역(0D, 초) → 유체(분포·시스, EEDF 가정) → 하이브리드(EM + 유체 + 운동론 전자) → PIC/MCC(EEDF·IEDF, 비쌈). 크누센 수 $\lambda/L$ 가 선택의 첫 잣대다.
표류-확산 모델: $\partial_tn + \nabla\cdot\Gamma = S$ , $\Gamma = \pm\mu nE - D\nabla n$ , 포아송과 결합. 반암시적 포아송은 유전 이완 시간 제한을, 샤페터-거멜 $\Gamma = (D/h)[B(-a)n_i - B(a)n_{i+1}]$ 은 큰 페클레 수의 진동을 없앤다.
LFA는 반응을 국소 $E/p$ 로, LMEA는 평균 에너지 방정식의 $\bar\varepsilon$ 로 정한다. 에너지 이완 길이가 장 변화 길이보다 길면(저압, 시스-글로 경계) LFA가 원천 위치를 틀린다.
PIC/MCC 사이클: 할당 → 장 풀이 → 보간 → 밀기 → 충돌(널 충돌법). 조건: $\Delta x \lesssim \lambda_D$ (아니면 유한 격자 가열), $\omega_{pe}\Delta t \lesssim 0.2$ (2를 넘으면 발산), $v\Delta t \lt \Delta x$ , 셀당 입자 수가 잡음을 정함.
1D PIC로 플라즈마 진동(보옴-그로스, 립프로그 분산), 이류 불안정(γ ≈ 0.35 ω_pe), 충돌 감쇠, 수치 가열을 직접 볼 수 있다.
ICP 하이브리드는 전자기(표피 깊이) ↔ 전자 에너지(이온화 원천) ↔ 유체 수송(밀도)을 수렴할 때까지 반복한다. 코일 분포가 웨이퍼 플럭스 균일도를 정한다.
교정과 검증을 나누고, 진단 오차를 비교에 넣는다. 진단의 계통 오차는 모델 매개변수에 흡수되어 다른 조건 예측을 틀리게 한다.
PIC 비용 ∝ 격자 × 입자 × 단계, 3D에서 대략 $n^2$ 로 늘어난다. 기준 ICP를 명시적 2D PIC로 푸는 것은 비현실적이라 유체·하이브리드가 장비 설계의 주력이다.

확인 퀴즈

1. 3 Torr PECVD 평행판(간격 10 mm)과 5 mTorr ICP 중, 국소장 근사 유체 모델이 더 잘 맞는 쪽과 이유는?

LFA는

\lambda_\varepsilon \ll L

일 때 맞다. 수 Torr에서

\lambda_\varepsilon

은 수십 µm–mm, 수 mTorr에서는 수 cm–챔버 크기로 EEDF가 비국소적이다.

2. 시스 안 격자 간격 0.25 mm에서 전자 표류 속도 10⁶ m/s, 확산 계수 50 m²/s이다. 셀 페클레 수와 적절한 이산화는?

a = vh/D = 10^6\times2.5\times10^{-4}/50 = 5 \gt 2

. 중앙 차분의 정상 해는

|a| \gt 2

에서 진동(음의 밀도)하고, 이는 시간 간격과 무관한 공간 이산화 문제다.

3. 1D 표류-확산 모델에서 포아송을 이전 단계 밀도로(명시적으로) 풀었더니 Δt = 1 ns에서 발산했다. 원인은?

명시적 결합에서 공간 전하는

e^{-t/\tau_d}

로 사라져야 하는데 Δt ≫ τ_d면 부호가 뒤집히며 커진다. 다음 단계 표류 전류를 포아송에 넣는 반암시적 방법으로 해결한다.

4. nₑ = 10¹⁶ m⁻³, Tₑ = 2 eV 플라즈마를 1D 명시적 PIC로 풀 때 권장 격자 간격과 시간 간격에 가장 가까운 것은?

\lambda_D = 7430\sqrt{2/10^{16}} \approx 1.05\times10^{-4}

\omega_{pe} = 56.4\sqrt{10^{16}} \approx 5.6\times10^9

s⁻¹.

\Delta t = 0.2/\omega_{pe} \approx 36

ps.

5. 열 플라즈마 PIC에서 계산을 빨리 하려고 Δx를 λ_D의 4배로 키웠더니 시간이 갈수록 Tₑ가 올랐다. 무엇이 일어났나?

\Delta x \gg \lambda_D

이면 격자보다 짧은 파장의 전하 요동이 긴 파장 모드와 섞여 에너지가 입자로 흘러든다. λ_D가 Δx/π 정도까지 커지면 가열이 멈춘다.

6. ICP 하이브리드 모델에서 전력을 2배로 올렸더니 흡수가 창 바로 아래로 더 몰렸다. 모델 안에서 일어난 고리는?

무충돌 표피 깊이

\delta = c/\omega_{pe} \propto n^{-1/2}

. 4×10¹¹ cm⁻³에서 약 8 mm다. 전자기 모듈과 유체 모듈이 밀도를 주고받으며 이 비선형 고리를 수렴시킨다.

7. 탐침 지름 분포에 모델을 교정했더니 χ²가 아주 작았다. 그런데 다른 압력 조건 예측이 25 % 틀렸다. 가장 그럴듯한 원인은?

교정은 측정과 모델의 차이를 매개변수로 메운다. 측정 자체에 계통 오차가 있으면 매개변수가 틀린 값으로 간다. 교정에 쓰지 않은 조건과 원리가 다른 진단으로 검증해야 드러난다.