Statistical Analysis of Electron-Defect Interaction in Thermoelectric Materials for the Development of Electron-Defect Averaged (EDA) Approach

Abstract

최근 인간 활동 등으로 인한 온실가스 배출량이 급증함에 따라 지구 온난화를 수반하는 기후변화 문제가 심각해지고 있다. 이러한 기후변화 문제를 해결하기 위해서는 화석연료 발전을 이용한 에너지 생산에서 벗어나 지속가능한 에너지 기술을 개발하고 기존 에너지 시스템의 효율을 증가시키는 것이 필수적이다. 에너지 효율을 증가시키는 저탄소 친환경 기술 가운데 각광받는 기술로는 열전 에너지 변환기술이 있다. 열전기술은 폐열을 전기로 변환시켜 에너지를 효율적으로 이용할 수 있도록 하는 기술이지만 높은 설치 단가와 낮은 발전 효율 때문에 널리 상용화되지 못하고 있다는 문제점을 가지고 있다. 따라서 이를 해결하기 위해 다양한 열전물질을 연구하고 조건을 최적화해 열전발전기의 성능을 높이는 연구가 필수적이다. 본 연구에서는 열전물질 가운데 비스무트 구리 옥시설파이드(bismuth copper oxysulfide: BiCuOS) 내에 생기는 구리 빈자리 결함(copper vacancy)에 의한 전자 산란현상을 측정하기 위해 Electron-Defect Averaged (EDA) 근사법을 개발하고 이를 검증하고자 하였다. BiCuOS는 결정질(crystalline)물질임에도 불구하고 구리의 약한 결합 때문에 낮은 열전도도(thermal conductivity)를 보인다는 장점이 있어 최근 열전물질로 각광받고 있다. 특히 BiCuOS에는 구리 결함이 일정한 농도로 생긴다는 사실이 알려져 있고 구리 결함으로 인해 생긴 정공(hole)이 전하운반체(charge carrier) 역할을 하며 p-type 반도체로 사용된다. 또한 구리 결함은 물질 내에서 전자를 산란시키거나 붙잡으면서 전자 이동도(carrier mobility)에 많은 영향을 주기 때문에 전자와 결함 간 상호작용(electron-defect interaction)을 계산하는 것이 열전물질의 성능평가에 매우 중요하다. EDA 근사법은 운동량에 대한 식인 M값을 에너지에 대한 식으로 바꿈으로써 BiCuOS와 같이 복잡한 물질에 대해서도 효과적으로 결정 결함에 의한 산란현상을 계산할 수 있도록 한 계산방법이다. M값은 결정 결함에 의한 전자 산란현상이 얼마나 많이 일어나는가를 나타내는 지표로 M^2값이 클수록 전자 산란이 많이 일어났음을 의미한다. EDA 근사법의 유효성을 검증하기 위해 본 연구에서는 6가지 크기의 슈퍼셀(supercell)을 모델링해 M^2값을 계산하였고 일정한 농도의 결정 결함 하에서 슈퍼셀의 크기에 관계없이 M^2값이 일정하게 계산될 것이라는 가설을 세웠다. BiCuOS의 전자구조를 분석하기 위해 제1원리계산(first principles calculation)의 일종인 밀도범함수이론(density functional theory: DFT)을 수행하였고 구리 결정 결함에 의한 전자 수송 특성(electronic transport property)의 변화를 계산하기 위해 볼츠만 수송 방정식(Boltzmann transport equation)을 이용해 제벡계수(Seebeck coefficient)와 전기 전도도(electrical conductivity)를 계산하였다. 첫번째로, 구리 결함을 포함하지 않은 결정구조(pristine unit cell)와 구리 결함을 포함한 결정구조(copper-deficient unit cell)의 밴드구조(band structure)와 상태밀도(density of state)를 계산한 결과, 구리 결함을 포함한 결정구조가 2배 이상 큰 띠 간격(bandgap)을 가지는 것으로 나타났다. 이는 구리가 가전자대에서 황과 반결합상태(antibonding state)를 이루면서 높은 에너지 상태를 점유하기 때문에 구리 결함이 생겼을 때 띠 간격이 넓어진 것이다. 또한 BiCuOS는 높은 제벡계수를 의미하는 heavy band와 높은 전기전도도를 나타내는 light band가 조합되어 열전물질로서 유리한 조건을 가짐을 알 수 있다. 두번째로, 가설 검증을 위해 가전자대(valence band)에서 계산된 M^2분포의 기술통계 분석 결과, 연구 초반에 세운 가설과 달리 M^2값의 평균이 슈퍼셀 크기의 1.7제곱에 반비례해서 작아진다는 결과를 얻었다. 이러한 결과는 첫째, M^2값이 슈퍼셀 사이즈가 충분히 커졌을 때 특정 값에 수렴할 가능성을 시사하고 둘째, EDA 근사법이 결정 결함에 의한 전자 산란현상을 제대로 계산하지 못했을 가능성을 의미한다. 그 원인으로는 BiCuOS가 이방성(anisotropic)구조를 갖는 물질이기 때문에 전자 산란과정에서 운동량의 중요성이 커지는데 EDA 근사법은 운동량을 무시했기 때문에 M^2값이 예상보다 작게 계산되었을 것이라고 예상할 수 있다. 따라서 추후 연구에서는 충분히 큰 슈퍼셀을 모델링해 M^2값의 수렴여부를 파악하거나 등방성(isotropic)구조를 갖는 물질에 대해 M^2값을 계산함으로써 EDA 근사법의 유효성을 검증할 수 있을 것이다. 세번째로, 정규성 검정을 통해 가전자대와 전도대(conduction band)의 M^2분포를 분석한 결과 모든 슈퍼셀에서 M^2값이 정규분포를 따르지 않는다는 결과를 얻었다. 특히 전도대의 M^2값은 극심한 왜도(skewness)와 이상치(outlier)를 보였고 정규성을 높이기 위해 이상치 제거(outlier removal)를 수행하였으나 정규성 검정 결과 여전히 정규분포를 따르지 않는다는 결과를 얻었다. 이는 첫째, 가전자대에 비해 전도대의 파동함수(wave function)이 더 복잡하기 때문에 본 연구에서 수행한 DFT 계산의 해상도(resolution)로는 이를 제대로 계산하지 못했고 둘째, 같은 에너지 범위 내에서 전도대의 DOS가 더 작기 때문에 전자와 구리 결함 간 산란을 일으킬 수 있는 에너지 상태의 수가 적어서 계산된 M^2값의 샘플 사이즈가 작아 대표성이 떨어지기 때문이다. 따라서 추후 연구에서 DFT 계산의 성능을 높여 복잡한 파동함수를 제대로 고려하고 더 넓은 에너지 범위에서 샘플링을 진행해 대표성을 띠는 M^2값을 계산할 필요가 있다. 다음으로, M^2분포의 정규성을 높이기 위해 거듭제곱 지수(power exponent)를 곱해 transformation을 시행한 결과, 가전자대의 2×2×1 슈퍼셀과 전도대의 2×2×1, 3×3×1 슈퍼셀의 분포가 정규분포를 따른다는 결과를 얻었다. 이를 통해 M^2분포가 z축 방향으로의 슈퍼셀 크기 증가에 x축과 y축으로의 변화보다 더 민감하게 변화함을 알 수 있다. 이를 바탕으로 가우시안 과정 회귀 모델(Gaussian Process Regression)을 이용해 M^2를 에너지에 대한 함수로 피팅해본 결과, transformation을 수행한 M^2분포가 가장 나은 피팅을 보여주었다. 그러나 이러한 결과는 M^2값의 샘플 크기가 작다는 점을 고려했을 때 신뢰성이 떨어진다. 따라서 추후 연구에서는 보다 많은 k-point에서의 샘플링을 수행하여 대표성이 있는 M^2샘플을 구하는 것이 중요하다. 마지막으로, 구리 결함에 의한 전자 수송 특성(electronic transport property)의 변화를 측정하고 EDA 근사법의 유효성을 검증하기 위해 일정 완화시간 근사법(constant relaxation time approximation: τ_constant)과 M^2값을 통해 직접 도출한 결정 결함으로 인한 완화시간 함수(relaxation time function: τ_constant)의 결과를 비교하였다. 그 결과, 제벡계수는 완화시간 근사법의 종류에 큰 영향을 받지 않았고 구리 결함의 농도와 음의 상관관계(negative correlation)를 보였다. 반면 전기전도도의 경우 완화시간 근사법의 종류에 따라 다른 양상을 보였는데 τ_constant를 적용한 결과, 전기전도도와 구리 결함 농도가 양의 상관관계(positive correlation)를 보였다. 반면 EDA 근사법을 통해 본 연구에서 도출한 τ_vacancy함수를 적용한 결과 전기전도도와 구리 결함 농도는 음의 상관관계를 보였는데 이는 전하 운반체의 증가에도 불구하고 τ_vacancy의 값이 이를 상쇄할 정도로 크게 감소하기 때문이다. 본 연구에서 M^2값이 부정확하게 계산됨에 따라 τ_vacancy는 τ_total에 기여하는 바가 매우 작게 계산되었기 때문에 EDA 근사법을 통해 얻은 τ_vacancy를 통해서는 물질의 전자 수송 특성을 계산하는 것이 어렵다는 결론을 얻었다. 이후, 적절한 M^2값이 계산되었다는 가정 하에 M^2값에 10^3을 곱해 다시 τ_vacancy를 계산하였고 그 결과, τ_vacancy의 영향이 τ_total의 결과를 결정할 정도로 커짐을 알 수 있었다. 이러한 결과는 정확한 M^2값이 계산될 경우, EDA 근사법을 통해 물질의 전자 수송 특성을 적절히 계산해낼 수 있다는 가능성을 시사한다. 따라서 추후 연구에서 앞서 언급한 보완점들을 통해 보다 정확한 M^2값과 τ_vacancy를 계산함으로써 EDA 근사법의 유효성을 검증할 수 있을 것이다.;In this study, we have developed and validated Electron-Defect Averaged (EDA) approach to evaluate the electron scattering phenomenon induced by the copper vacancy in BiCuOS. EDA method efficiently evaluates the electron-defect scattering phenomena by converting the momentum-dependent electron-defect interaction matrix (M) into the energy-dependent function. The M^2 value indicates the degree to which electron scattering occurs. We hypothesize that the M^2 values would be independent of the supercell size under the constant defect concentration. To prove the hypothesis, we calculate M^2 values for six cases of supercell structure and investigate the M^2 distributions using diverse statistical analysis tools. Furthermore, the changed electrical transport properties induced by Cu vacancy are calculated to assess the validity of the EDA method. We employ DFT calculation to investigate the electronic structure of BiCuOS using the Quantum ESPRESSO software. Also, we calculate the Seebeck coefficient and electrical conductivity using the BoltzTraP code which is based on the Boltzmann Transport Theory. The constant relaxation time and relaxation time functions have been employed to calculate the transport properties. First, the Cu-deficient unit cell shows a larger bandgap value than that of the pristine unit cell. This is caused by the Cu-S antibonding states which lead to high energy states near the valence band maximum (VBM) and thus, Cu vacancy widens the bandgap. Also, the combination of heavy and light bands near the VBM indicates desirable thermoelectric performance. The heavy band denotes a large Seebeck coefficient which has been caused by sulfur and the light band shows high electrical conductivity induced by Cu vacancies. Next, we analyze the descriptive statistics of valence band M^2 distributions. As a result, the M^2 mean values are inversely proportional to the 1.7 square of the cell size which goes against the initial hypothesis. It indicates the possibility that M^2 may converge to a certain value in the larger supercell. Therefore, to check the convergence of M^2, large enough supercells should be modeled for the M^2 calculation. Moreover, another reason for the M^2 values different from the hypothesis is the anisotropy of BiCuOS structures. For anisotropic materials, the electron momentum can play a significant role in the electron-defect scattering phenomena. Thus, neglecting momentum can deter the EDA method from fully estimating the scattering phenomena and lead to inaccurate M^2 values. Then, we conduct a normality test of M^2 distributions and the results show that M^2 datasets from all the supercell structures are not normally distributed. Especially, the extreme asymmetry of CB M^2 distribution is mainly caused by fewer scattering events in CB than in VB. To enhance the normality, we remove outliers and conduct transformation. Although the outlier removal does not change the overall trend, the transformation significantly increases the normality of M^2 distribution. The transformed M^2 distributions of 2×2×1 structure VB and CB show normality and the 3×3×1 structure CB M^2 distribution after the transformation and outlier removal also shows normality according to the Shapiro-Wilk test. However, it is hard to rely on the normality test results because of the small sample size of CB M^2 datasets. After that, we conduct Gaussian Process Regression (GPR) to fit the M^2 datasets and obtain the predictive mean function more easily. GPR for transformed VB M^2 distributions presents better fitting with shorter length scales. However, the overall fitting seems not successful given the M^2 values that are widely spread out from the mean function. To obtain a better fitting, it is necessary to increase the number of sampling k-points for the M^2 calculation. Finally, we calculate the Seebeck coefficient and electrical conductivity using the τ_constant, τ_vacancy and τ_total. τ_vacancy is calculated from M^2 values and described as a regressed function of energy. The Seebeck coefficients show similar results from each of the τ, indicating they are rarely affected by τ values. The Seebeck coefficients have a negative correlation with the doping level and a positive correlation with temperature. The electrical conductivity calculated from τ_constant results in a positive correlation with the doping level while τ_vacancy brings about a negative correlation with the doping level. However, since τ_vacancy rarely contributes to the results of τ_total, its impacts could be neglected. This result indicates that it is hard to evaluate the proper transport properties using the τ_vacancy from our inaccurate M^2 results. Therefore, we recalculate the transport properties assuming we obtain reasonable M^2 values from EDA. As a result, the modified τ_vacancy dominates τ_total, indicating that if the EDA method calculates proper M^2 values, we can somehow evaluate the transport properties using the τ_vacancy only. Therefore, it is necessary to obtain accurate τ_constant and proper M^2 values to verify the EDA approach and assess the transport properties in further study.