Estudio de propiedades físicas del silicio mediante el uso de la termodinámica estadística
Full text
Equation Chapter 1 Section 1 Trabajo Fin de Máster Máster Universitario en Diseño Avanzado en Ingeniería Mecánica Estudio de propiedades físicas del silicio mediante el uso de la termodinámica estadística Autor: Carlos Sánchez Martín Tutora: Pilar Ariza Moreno Dep. Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016
Trabajo Fin de Máster Máster Universitario en Diseño Avanzado en Ingeniería Mecánica Estudio de propiedades físicas del silicio mediante el uso de la termodinámica estadística Autor: Carlos Sánchez Martín Tutora: Pilar Ariza Moreno Catedrática de Universidad Dep. de Mecánica de Medios Continuos y Teoría de Estructuras Escuela Técnica Superior de Ingeniería Universidad de Sevilla Sevilla, 2016
Trabajo Fin de Máster: Estudio de propiedades físicas del silicio mediante el uso de la termodinámica estadística Autor: Carlos Sánchez Martín Tutora: Pilar Ariza Moreno El tribunal nombrado para juzgar el Proyecto arriba indicado, compuesto por los siguientes miembros: Presidente: Vocales: Secretario: Acuerdan otorgarle la calificación de: Sevilla, 2016 El Secretario del Tribunal
iii Agradecimientos A T´ ı, que siempre es´ as ah´ ı, a´ un cuando no me doy cuenta. Septiembre 2016, Sevilla (Espa˜ na)
0. ´ Indice general ´ Indice general iv 1. Introducci´ on 1 1.1. Simulaciones moleculares . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.2. Materiales termoel´ ectricos ............................ 2 2. El silicio 4 2.1. Descripci´ on de la red cristalina del silicio . . . . . . . . . . . . . . . . . . . 4 2.2. Conductividad t´ ermicadelsilicio ........................ 5 3. Teor´ ıa termodin´ amica estad´ ıstica de no equilibrio 7 3.1. Aproximaci´ on de m´ axima entrop´ ıa y relaciones de equilibrio . . . . . . . . 7 3.2. Relaciones cin´ eticas................................ 9 3.2.1. Balance local y desigualdad de disipaci´ onlocal............ 9 3.2.2. Cin´ eticadeOnsager............................ 10 3.2.3. Cin´ eticaLineal .............................. 11 3.3. Teor´ ıa variacional de campo medio . . . . . . . . . . . . . . . . . . . . . . . 11 3.4. Evaluaci´ on de las medias de fase de campo medio . . . . . . . . . . . . . . . 13 3.4.1. Expansi´ onmultipolos........................... 13 3.4.2. Cuadratura de Gauss . . . . . . . . . . . . . . . . . . . . . . . . . . . 14 iv
´ INDICE GENERAL v 4. Modelo simple para el estudio de la conductividad t´ ermica en nanohilos de silicio 16 5. M´ etodos num´ ericos para la optimizaci´ on de problemas no lineales 21 5.1. Gradiente conjugado no lineal . . . . . . . . . . . . . . . . . . . . . . . . . . 22 5.2. Relajaci´ on din´ amicaadaptativa ......................... 24 6. Potenciales termodin´ amicos para el silicio 28 6.1. Potencial Stillinger-Weber . . . . . . . . . . . . . . . . . . . . . . . . . . . . 28 6.2. PotencialMEAM.................................. 30 6.3. Comparaci´ on entre la aproximaci´ on de multipolos y la cuadratura de Gauss para el potencial Stillinger-Weber . . . . . . . . . . . . . . . . . . . . . . . . 32 6.4. Comparaci´ on entre potenciales en t´ erminos energ´ eticos ........... 35 6.5. Validaci´ on de los potenciales para representar el comportamiento del silicio 40 6.5.1. Matriz de constantes de fuerza . . . . . . . . . . . . . . . . . . . . . 40 6.5.2. Matriz din´ amica.............................. 43 6.5.3. Dispersi´ ondefonones .......................... 44 7. Conclusiones y trabajos futuros 49 A. Ap´ endice A 50 A.1. Derivadas del potencial Stillinger-Weber . . . . . . . . . . . . . . . . . . . . 50 A.1.1. Derivadas del t´ ermino vij(rij) ...................... 50 A.1.2. Derivadas del t´ ermino vij(rji) ...................... 51 A.1.3. Derivadas del t´ ermino vjik(rij,rik).................... 51 A.1.4. Derivadas del t´ ermino vjik(rji,rjk).................... 51 A.2. Expresi´ on del t´ ermino ∆jVi(q)paraelSW ................... 52 B. Ap´ endice B 54 B.1. Derivadas del potencial MEAM . . . . . . . . . . . . . . . . . . . . . . . . . 54 B.1.1. Derivadas del t´ ermino Sij ........................ 55 B.1.2. Derivadas del t´ ermino ρ(0) i........................ 56 B.1.3. Derivadas del t´ ermino ρ(k) ij ........................ 57
´ INDICE GENERAL vi B.1.4. Derivadas del t´ ermino t(k) i........................ 58 B.1.5. Derivadas del t´ ermino Γ(k) i........................ 59 B.1.6. Derivadas del t´ ermino ρi......................... 59 B.2. Expresi´ on del t´ ermino ∆jVi(q)paraelMEAM ................. 60 Bibliografa i ´ Indice de figuras i Bibliograf´ ıa iii
1. Introducci ´ on 1.1. Simulaciones moleculares El objetivo de este trabajo es obtener propiedades mec´ anicas y t´ ermicas mediante t´ ecnicas computacionales de nanoestructuras como los nanohilos de silicio. A diferencia del mundo macrosc´ opico, cuando se quieren conocer propiedades de materiales a escala microsc´ opica su determinaci´ on resulta mucho m´ as compleja al abandonar el mundo continuo y pasar a trabajar con sistemas discretos, adem´ as de que las propiedades que se miden durante la simulaci´ on no son directamente las que se quieren conocer y por tanto hay que establecer relaciones entre unas y otras que sean consistentes con las leyes f´ ısicas. La mec´ anica estad´ ıstica [1] proporciona la base te´ orica para estudiar el comportamiento termomec´ anico de la muestra de estudio. A trav´ es de un tratamiento probabil´ ıstico, es posible conocer observables macrosc´ opicos de la muestra a partir del conocimiento de sus variables microsc´ opicas. Una de las metodolog´ ıas m´ as usadas para calcular estos observables es el denominado m´ etodo de Monte Carlo [2], y consiste en evaluar de forma aleatoria posibles estados del sistema, de modo que cuantos m´ as estados aleatorios tengamos y cuanto m´ as se acerquen estos al estado real en el que se encuentra el sistema, menor ser´ a el error cometido. Otra de las metodolog´ ıas para realizar simulaciones moleculares m´ as extendida es la din´ amica molecular [3]. Esta t´ ecnica consiste en resolver las ecuaciones de movimiento de Newton para un sistema de Npart´ ıculas hasta que las propiedades del sistema no var´ ıen, momento en el que se habr´ a alcanzado el equilibrio. Una vez alcanzado el equilibrio, se puede medir cualquier observable de la muestra, como la temperatura o la energ´ ıa, a partir del conocimiento de la posici´ on y velocidad de cada part´ ıcula. En este trabajo se utilizar´ a una teor´ ıa termodin´ amica estad´ ıstica de no equilibrio [4], que permitir´ a el estudio de sistemas fuera del equilibrio, pudiendo trabajar en todo momento con valores macrosc´ opicos de las variables.
1. Introducci ´ on 2 1.2. Materiales termoel´ ectricos Aproximadamente el 90% de la producci´ on de energ´ ıa del mundo es generada por motores t´ ermicos que usan combustibles f´ osiles como fuente de calor y que operan generalmente entre un 30 −40% de eficiencia, de modo que cerca de 15 teravatios de calor se pierde al ambiente. Los materiales termoel´ ectricos pueden convertir energ´ ıa t´ ermica en energ´ ıa el´ ectrica o viceversa y podr´ ıan convertir parte de este calor perdido a electricidad aprovechable. La eficiencia de un material termoel´ ectrico se caracteriza por su figura de m´ erito (ZT), un par´ ametro adimensional que viene dado por ZT =S2GT κ(1.1) donde Ses el coeficiente de Seebeck, Ges la conductancia electr´ onica, Tes la temperatura y κes la conductividad t´ ermica, respectivamente. Los materiales que tienen ZT ≈1 se consideran como buenos termoel´ ectricos, mientras que ZT > 3 se requiere para competir con refrigeradores convencionales o generadores. Un material termoel´ ectrico ideal es aquel que se comporta con los fonones como vidrio y los dispersa, dando un valor bajo de κ, mientras que con los electrones es como un cristal perfecto sin dispersi´ on, dando un valor alto del factor de potencia (S2G). El material termoel´ ectrico m´ as ampliamente usado comercialmente es el telururo de bismuto (Bi2T e3) y sus aleaciones con Sb,Se, etc, que tienen una figura de m´ erito ZT ≈1. Es dif´ ıcil escalar el telururo de bismuto macizo para conversi´ on de energ´ ıa a gran escala, pero fabricar nanoestructuras sint´ eticas con este prop´ osito es incluso m´ as dif´ ıcil y caro. El silicio, por otro lado, es el conductor m´ as abundante y usado, con una gran infraestructura industrial para su procesamiento a bajo coste. El silicio macizo, sin embargo, tiene una figura de m´ erito ZT ≈0,01 a 300K[5]. La distribuci´ on espectral de fonones que contribuyen a la κdel silicio a temperatura ambiente es bastante amplio. Debido a que el ratio de dispersi´ on fon´ on-fon´ on Umklapp se escala como ω2, donde ωes la frecuencia del fon´ on, los fonones ac´ usticos de baja frecuencia (o gran longitud de onda) tienen un gran camino libre medio y contribuyen significativamente a la κa altas temperaturas [6, 7] . As´ ı, incorporando elementos de dispersi´ on de fonones a distintas escalas se puede esperar una alta disminuci´ on de la κdel silicio. Los nanohilos de silicio son candidatos prometedores como materiales termoel´ ectricos eficientes [8, 9, 10] a pesar de que el silicio macizo no lo sea. Los nanohilos conducen bien la carga, pero la conductividad t´ ermica se reduce dr´ asticamente debido a la disper- si´ on de fonones en las superficies rugosas. Este trabajo se organiza de la siguiente manera: en el Cap´ ıtulo 2 se hace una descripci´ on general del material a estudiar, el silicio. En el Cap´ ıtulo 3 se presenta la teor´ ıa termodin´ amica estad´ ıstica de no equilibrio que se utilizar´ a para el estudio. Seguidamente, en el Cap´ ıtulo 4 se presenta un modelo simple de c´ alculo de conductividad t´ ermica para nanohilos de silicio. En el Cap´ ıtulo 5 se comentar´ an los m´ etodos matem´ aticos utilizados para la resoluci´ on de los sistemas de ecuaciones no lineales planteados. En el Cap´ ıtulo 6 se introducir´ an y describir´ an los potenciales interat´ omicos utilizados para simular las fuerzas entre ´ atomos y generar un modelo m´ as complejo que incluya la parte mec´ anica.
3. Teor´ ıa termodin ´ amica estad´ ıstica de no equilibrio 9 gran can´ onica de Gibbs fuera del equilibrio. Adem´ as, por la misma analog´ ıa se puede interpretar θi=1 kBβi (3.9) como la temperatura absoluta de la part´ ıcula i, que puede variar de una part´ ıcula a otra. Adem´ as, la entrop´ ıa total del sistema se puede identificar con la entrop´ ıa de la informaci´ on para la distribuci´ on gran can´ onica S=kB{β}T{e}+kBlogΞ(3.10) De forma directa se puede calcular la relaci´ on local de equilibrio βi=1 kB ∂S ∂ei ({e}) (3.11) que relaciona la energ´ ıa de cada part´ ıcula con su temperatura. 3.2. Relaciones cin´ eticas Una vez conocidas las expresiones que permiten conocer los valores macrosc´ opicos de las variables del sistema considerado, se procede a introducir la formulaci´ on que se utilizar´ a para el estudio de procesos cin´ eticos. 3.2.1. Balance local y desigualdad de disipaci´ on local Se comienza examinando el balance de energ´ ıa a nivel de part´ ıcula. La energ´ ıa de la part´ ıcula ipuede identificarse como el valor esperado del hamiltoniano de la part´ ıcula ui=hhii(3.12) para i= 1,...,N. Sup´ ongase, adem´ as, que los hamiltonianos locales hidependen de variables macrosc´ opicas. Entonces, el balance de energ´ ıa a nivel de part´ ıcula se escribe como ˙ ui=˙ wi+ri(3.13) donde ˙ wies el trabajo mec´ anico externo y ries el flujo de calor en la part´ ıcula i. Suponiendo que el flujo de calor rien la part´ ıcula ies de la forma ri=X j,i Rij (3.14) donde la suma se extiende a todas la part´ ıculas distintas de iy
3. Teor´ ıa termodin ´ amica estad´ ıstica de no equilibrio 10 Rij =−Rji (3.15) es el flujo de calor discreto desde la part´ ıcula jhasta la i. Ahora, se busca formular el balance de entrop´ ıa y una desigualdad de disipaci´ on a nivel de part´ ıcula. Para ello, se considera cada part´ ıcula como un sistema en equilibrio con su entorno, considerado como un ba˜ no de calor, y se espera que aparezca disipaci´ on como resultado del flujo de calor entre part´ ıculas. Como analog´ ıa con la termodin´ amica continua, se puede definir la variaci´ on de entrop´ ıa a nivel de part´ ıcula por medio de la relaci´ on ˙ si kBβi =˙ ui−˙ wi(3.16) Adicionalmente, si se considera un subsistema formado por dos part´ ıculas, iyj, la cantidad Σij =1 θi ˙ ui−˙ wi−X k,i,j Rik +1 θj ˙ uj−˙ wj−X l,i,j Rjl =Rij θi +Rji θj =PijRij (3.17) donde se escribe Pij =1 θi−1 θj =kB(βi−βj) (3.18) mide la tasa de producci´ on interna de entrop´ ıa para el par de part´ ıculas, es decir, la tasa de entrop´ ıa dentro del par en exceso de la tasa de entrop´ ıa debido al calor externo suministrado. Entonces, se postula la desigualdad de producci´ on interna local de entrop´ ıa Σij ≥0 (3.19) La Ecuaci´ on 3.19 es an´ aloga a la desigualdad de disipaci´ on cl´ asica de Clausius- Duhem de la termodin´ amica continua. 3.2.2. Cin´ etica de Onsager Evidentemente, la desigualdad local de disipaci´ on (3.19) restringe el tipo de proceso, o de evoluci´ on, de la temperatura que puede ocurrir. Adem´ as, la desigualdad (3.19) establece el emparejamiento dual entre flujos y fuerzas: los flujos de calor Rij est´ an impulsados por el gradiente dual discreto de temperatura Pij. As´ ı, por analog´ ıa con la teor´ ıa
3. Teor´ ıa termodin ´ amica estad´ ıstica de no equilibrio 11 de Onsager de relaciones cin´ eticas, la desigualdad local de disipaci´ on (Ecuaci´ on 3.19) sugiere una ley cin´ etica del tipo Rij =∂Ψ ∂Pij ({P}) (3.20) donde Ψ({Pij}) es un potencial cin´ etico discreto. La Ecuaci´ on 3.20 se puede considerar como una ley discreta de Fourier de conducci´ on de calor. La ecuaci´ on discreta de conservaci´ on de energ´ ıa puede escribirse entonces como ˙ si kBβi =X j,i ∂Ψ ∂Pij ({P}) (3.21) 3.2.3. Cin´ etica Lineal La relaci´ on cin´ etica considerada es una relaci´ on emp´ ırica, que tiene que ser modelada de la misma forma que los potenciales interat´ omicos. Para obtener un modelo cin´ etico simple con el que trabajar, se considerar´ a que la relaci´ on entre flujos y fuerzas es lineal y que el potencial cin´ etico es la suma de las contribuciones de cada enlace entre part´ ıculas. De esta forma, el potencial tiene la forma general Ψ({P}) = X hi,ji Aijθ2 ijP2 ij (3.22) donde hi,jirepresenta el enlace definido por las part´ ıculas iyjy la suma se extiende sobre todos los enlaces del sistema. Adem´ as, se ha introducido la temperatura media de enlace definida como θij =1 2(θi+θj) (3.23) El coeficiente Aij es un coeficiente cin´ etico que depende, de forma general, de la temperatura y es diferente para cada material. Entonces, la ley discreta de Fourier ( Ecuaci´ on 3.20) se escribe como ˙ si kBβi =X hi,ji Aijθ2 ijPij =X hi,ji Aijθ2 ijkB(βi−βj) (3.24) 3.3. Teor´ ıa variacional de campo medio A pesar de la simplicidad formal de la teor´ ıa mec´ anica estad´ ıstica del no equilibrio presentada, el c´ alculo de potenciales termodin´ amicos en forma cerrada es, en general, intratable. Para lidiar con este problema, consid´ erese un conjunto de funciones de densidad
3. Teor´ ıa termodin ´ amica estad´ ıstica de no equilibrio 12 de probabilidad de prueba P0yp∗ 0la mejor probabilidad dentro del conjunto obtenida por medio del principio de m´ axima entrop´ ıa Lp∗ 0,{β}= sup p0∈P0L[p0,{β}](3.25) para {β}fijo. En la pr´ actica, se elige una funci´ on de densidad de prueba de la forma ρ0=1 Ξ0 e−{β}T{h0}(3.26) donde Ξ=1 h3NZΓ e−{β}T{h0}dqdp (3.27) y{h0}en alg´ un conjunto H0de hamiltonianos locales de prueba. Entonces, el principio de m´ axima entrop´ ıa (Ecuaci´ on 3.25) es equivalente a minimizar el funcional F[{h0},{β}]=kB{β}T{hh−h0i0}−kBlogΞ0(3.28) Se considerar´ a un hamiltoniano de prueba de la siguiente forma h0i({q},{p}) = 1 2mi|pi−pi|2+mω2 i 2|qi−qi|2(3.29) de modo que los hamiltonianos estar´ an parametrizados por los campos {q}= (qi)N i=1, {p}= (pi)N i=1, y {ω}= (ωi)N i=1, que pueden ser considerados como la posici´ on y el momento medios de las part´ ıculas y su frecuencia de vibraci´ on efectiva respectivamente. El valor ´ optimo de esos campos se obtiene variacionalmente mediante la minimizaci´ on del funcional. Adem´ as, n´ otese los hamiltonianos de prueba se han elegido de forma que sean cuadr´ aticos y que solo dependan de la posici´ on y el momento de la part´ ıcula. De esta forma, el funcional puede ser escrito de la forma F[{h0},{β}]=kB N X i=1 [βihhii0−3 + 3log(~βiωi)](3.30) Bajo condiciones cuasiest´ aticas, los valores ´ optimos {q∗},{p∗}, y {ω∗}se obtienen minimizando Fbajo {β}constante. Las correspondientes ecuaciones de Euler-Lagrange se escriben ∂F ∂qi =∂ ∂qikB{β}T{hhi0}= 0 ∂F ∂pi =∂ ∂pikB{β}T{hhi0}= 0 ∂F ∂ωi =∂ ∂ωikB{β}T{hhi0}+3kB ωi = 0 (3.31)
3. Teor´ ıa termodin ´ amica estad´ ıstica de no equilibrio 13 3.4. Evaluaci ´ on de las medias de fase de campo medio Como se puede observar en la Ecuaci´ on 3.30, la evaluaci´ on del funcional requiere el c´ alculo de los hamiltonianos locales hhii0, es decir, la media de fase de los hamiltonianos locales con respecto a la funci´ on de densidad de probabilidad de campo medio. Si, por ejemplo, los hamiltonianios locales son de la forma hi({q},{p}) = 1 2mi|pi|2+Vi({q},{p}) (3.32) donde Vi({q},{p}) es un potencial interat´ omico multiespecie, entonces se tiene hhii0=3 2βi +1 2mi|pi|2+hVii0(3.33) En general, la evaluaci´ on exacta de este valor efectivo no es posible, a´ un para un hamiltoniano de prueba simple como el propuesto. Para calcular entonces dicho valor, se aproxima la integral utilizando dos m´ etodos distintos: expansi´ on multipolos y cuadratura de Gauss. 3.4.1. Expansi ´ on multipolos Sea µuna medida acotada sobre Rny sea f∈C(Rn). La acci´ on de la medida µsobre fse escribe µ(f) = Zf(x)dµ(x) (3.34) La aproximaci´ on por multipolos de grado kde µes la medida µktal que µk(f) = X |α|≤k cαDαf(0) (3.35) donde αes un ´ ındice en N3y los coeficientes {cα,|α|≤k}son elegidos para que la aproximaci´ on sea exacta para polinomios de grado ≤k, tal que µk(xα) = µ(xα),|α|≤k(3.36) Este requisito implica que cα=1 α!Zxαdµ(x) (3.37) Utilizando esta metodolog´ ıa, la expansi´ on multipolos de segundo orden de hVii0tie- ne la forma
3. Teor´ ıa termodin ´ amica estad´ ıstica de no equilibrio 14 hVii0=ZVi({q}) N Y j=1 (~βjωj)3e−βj1 2mj|pj−pj|2−βj mjω2 j 2|qj−qj|2 dqdp h3N =ZVi({q}) N Y j=1 βjmjω2 j 2π 3 2 e−βj mjω2 j 2|qj−qj|2 dq ≈Vi({q}) + N X j=1 1 2βjmjω2 j ∆jVi({q}) (3.38) donde ∆jVi({q}) = 3 X α=1 ∂2 qjα Vi({q}) (3.39) De este modo, se puede considerar que el potencial est´ a termalizado en tanto en cuanto las temperaturas (βi= 1/kBθi) aparecen involucradas en el c´ alculo del valor efectivo del mismo. 3.4.2. Cuadratura de Gauss Sup´ ongase que cada funci´ on Viinvolucra un n´ umero finito de ´ atomos vecinos. Entonces, la integral hViipuede ser calculada usando la regla de cuadratura Hermite-Gauss apropiada para la dimensi´ on del espacio. Para el caso general de un potencial de interacci´ on de ncuerpos φ(q1,...,qn) el valor esperado resulta hφ(q1,...,qn)i=1 h3nZΓ φ(q1,...,qn) n Y i=1 ρidqidpi =1 h3nZΓ φ(q1,...,qn) n Y i=1 (~βiωi)3 n Y i=1 e−βi 2mi|pi−pi|2−βimiω2 i 2|qi−qi|2dqidpi =Qn i=1(~βiωi)3 h3n Zn Y i=1 e−βi 2mi|pi−pi|2dpi "Zφ(q1,...,qn)e−βimiω2 i 2|qi−qi|2dqi# =Qn i=1(~βiωi)3 h3n n Y i=1 r2πmi βi!3 "Zφ(q1,...,qn)e−βimiω2 i 2|qi−qi|2dqi# =Qn i=1(~βiωi)3 h3n n Y i=1 r2πmi βi!3 n Y i=1 1 ωir2 βimi!3 "Zφ(q1,...,qn)e−|xi|−...−|xn|dxi···dxn# = 1 √π!3n"Z˜ φ(x1,...,xn)e−|xi|−...−|xn|dxi···dxn# (3.40)
3. Teor´ ıa termodin ´ amica estad´ ıstica de no equilibrio 15 donde la Ecuaci´ on 3.40 es de dimensi´ on 3ny se ha hecho el siguiente cambio de variables xi=rβimi 2ωi(qi−qi) ˜ φ(x1,...,xn) = φ(q1(x1),...,qn(xn)) (3.41) Para una regla de cuadratura que utiliza Mpuntos la integral anterior se reduce a hφ(q1,...,qn)i≈1 π3n 2M X k=1 ˜ φ(x1,...,xn)wk(3.42) donde wkes el peso de cuadratura correspondiente. Debido al tama˜ no del espacio de integraci´ on que requiere un potencial multicuerpo como el MEAM, se utilizar´ a una aproximaci´ on de 3er orden. Para el caso de un potencial que considere nvecinos para calcular la energ´ ıa de un ´ atomo, se tendr´ an puntos de cuadratura de dimensi´ on m= 3nhasta un total de M= 2m. Para un espacio de dimensi´ on mlos puntos de cuadratura propuestos por [20] son ±rm 2,0,...,0! ... 0,0,...,±rm 2! (3.43) y los pesos se escriben como wk=πm 2 2mk= 1,...,M (3.44)
4. Modelo simple para el estudio de la conductividad t´ ermica en nanohilos de silicio A continuaci´ on se pretende reproducir los resultados experimentales obtenidos en [16]. Para ello, se utilizar´ a la teor´ ıa expuesta en el Cap´ ıtulo 3, pero realizando varias simplificaciones. En primer lugar, se considera que el sistema es r´ ıgido y por tanto los ´ atomos est´ an fijos, de modo que solo es necesario resolver la ecuaci´ on discreta de Fourier (Ecuaci´ on 3.21). Adem´ as, se considerar´ a que se trabaja en estado estacionario. Se define una celda computacional Ωconsistente en una rebanada de expesor lformada por 20 (111) planos at´ omicos y se utiliza un sistema de coordenadas cartesiano cuyo eje zcoincide con el eje del nanohilo. Adem´ as, se establece la temperatura de un ´ atomo a temperatura ambiente, denominada T0, y un gradiente t´ ermico Tz=∆T L(4.1) donde ∆T= 5.0 Kes aproximadamente la m´ axima variaci´ on de temperatura observada en el experimento, y L= 3.0 µm es la longitud aproximada del nanohilo. En r´ egimen estacionario, el campo de temperatura {θ}se obtiene resolviendo el problema m´ ın {θ}X hi,ji 1 2Aijθ2 ijP2 ij (4.2) s.a. X i∈Γ0 θiσi−X i∈Γ20 θiσi=|Ω|Tz(4.3) s.a. θi0=T0(4.4) donde Γ0yΓ20 son el primero y el ´ ultimo plano at´ omico de la rebanada respectivamente, yσies el rea de la secci´ on correspondiente al ´ atomo i, que viene dada por
4. Modelo simple para el estudio de la conductividad t´ ermica en nanohilos de silicio 17 σi=πR2 N(4.5) siendo Rel radio del nanohilo. Haciendo uso del m´ etodo de los multiplicadores de Lagrange, el problema se puede escribir como X hi,ji 1 2Aijθ2 ijP ij2−Jz T2 0 X i∈Γ0 θiσi−X i∈Γ20 θiσi−|Ω|Tz −λ(θi0−T0)=0 (4.6) donde Jzyλson los multiplicadores de Lagrange introducidos, representando el primero de ellos el flujo medio de calor que atraviesa el nanohilo. De este modo, la conductividad t´ ermica efectiva se puede calcular como κ=Jz Tz (4.7) Figura 4.1. Modelo computacional para la conducci´ on de calor en nanohilos de silicio (a) D = 56 nm; τa= 1.7 nm (b) D = 115 nm; τa= 3.5 nm. ´ Atomos pertenecientes a la capa amorfa coloreados en rojo. Como se puede ver en [16], los nanohilos usados durante el experimento ten´ ıan una capa rugosa de un espesor aproximado t= 1 nm para un di´ ametro de 22 nm. Se han realizado muchos estudios [21, 22, 23, 24, 25] sobre la efecto de esta capa amorfa superficial en la conductividad t´ ermica de los nanohilos de silicio. Seg´ un se ha observado en [26], la relaci´ on entre el espesor de esta capa rugosa y el radio del nanohilo es aproximadamente t≈0,06 R. En este estudio se ha utilizado esta relaci´ on, y dicha capa amorfa se genera desplazando de forma aleatoria los ´ atomos desde su posici´ on en la red perfecta hasta un m´ aximo del 10% del par´ ametro de red. En concreto, la posiciones at´ omicas medias se obtienen como qi= q(0) i+pidi,si i ∈Ωa; q(0) i,en otro caso (4.8)
4. Modelo simple para el estudio de la conductividad t´ ermica en nanohilos de silicio 18 donde q(0) irepresenta la posici´ on media del ´ atomo ien la red perfecta del Si y Ωa= ∈Ω|rq(0)2 ix +q(0)2 iy > R −τa (4.9) identifica a los ´ atomos localizados dentro de la capa amorfa. El espesor de la capa amorfa τase ha calculado como τa= 0,06 Rcomo se ha comentado antes, mientras que piydison un m´ odulo y una direcci´ on que se generan aleatoriamente. El m´ odulo se genera mediante una distribuci´ on Gausiana Npmax 2,(pmax 8)2, con pmax =a/10, mientras que la direcci´ on se genera con una distribuci´ on uniforme en una esfera unidad S2=nx∈R3|kxk2= 1o. En la Figura 4.1 se representan los modelos computacionales para nanohilos de D= 56 nm yD= 115 nm, donde la capa amorfa se colorea en rojo. En cuanto al coeficiente cin´ etico de transporte, se supone que es de la forma Aij = Aa,si i ∈Ωao j ∈Ωa,kqi−qjk< rc Ac,si i,j ∈Ω\Ωa,kqi−qjk< rc 0,si kqi−qjk< rc (4.10) con Aa(T ,D) = α(D)T3+β(D)T2+γ(D)T(4.11) +δ(D) (4.12) Ac(T ,D) = Aa(T ,D) 0,005 (4.13) donde AcyAason los coeficientes de transporte para las fases cristalina y amorfa y rc es un radio de corte. Se establece el valor de dicho radio seg´ un rc= (r1+r2)/2, donde r1 yr2son los radios de la primera y segunda capa de vecinos en la red perfecta. De este modo, solo se considera la interacci´ on con los primeros vecinos en la red perfecta. Para fijar la relaci´ on entre AcyAase han comparado los resultados presentados para el silicio cristalino en [5] con los obtenidos para el silicio amorfo en [27], y se ha observado que para el rango de temperaturas de inter´ es dicha relaci´ on se encuentra en torno al valor escogido. As´ ı mismo, los coeficientes α(D),β(D),γ(D) y δ(D) dependen del di´ ametro del nanohilo. Al igual que para la dependencia de AcyAacon la temperatura, para los coeficientes α(D),β(D),γ(D) y δ(D) tambi´ en se han usado polinomios c´ ubicos. En concreto, tienen las expresiones: α(D) = −(1,526 ·10−16)·D3+ (3,032 ·10−13)·D2 −(1,27 ·10−10)·D+ (1,538 ·10−8)(4.14) β(D) = (9,616 ·10−14)·D3−(1,909 ·10−10)·D2 + (7,565 ·10−8)·D−(8,689 ·10−6)(4.15)
5. M´ etodos num´ ericos para la optimizaci´ on de problemas no lineales 25 diferencias centrales con estaci´ on intermedia. De este modo, las velocidades se definen en el punto medio del paso de tiempo y la aproximaci´ on para las derivadas temporales son ˙ qn+1 2 a=1 hqn+1 a−qn a(5.11) ¨ qn a=1 h˙ qn+1 2 a−˙ qn−1 2 a(5.12) donde hes el incremento de tiempo fijado. El vector de velocidad para computar la fuerza amortiguadora se puede estimar como ˙ qn a=1 2˙ qn+1 2 a−˙ qn−1 2 a(5.13) Introduciendo las Ecuaciones 5.12 y 5.13 en 5.10 y reagrupando t´ erminos en la Ecuaci´ on 5.11, se llega a que la velocidad y el desplazamiento en el siguiente paso de tiempo son ˙ qn+1 2 a= 2−ch 2 + ch!qn−1 2 a+ 2h 2 + ch!M−1Fn ext −Fn int qn+1 a=qn a+h˙ qn+1 2 a (5.14) donde la inversa de la matriz Mse calcula de forma trivial al ser Muna matriz diagonal. La ecuaci´ on anterior nos da el proceso iterativo para calcular qa, pero el paso inicial no est´ a determinado. Dado que ˙ q0 ayq0 ason conocidos, se puede obtener la siguiente relaci´ on ˙ q 1 2 a=h 2M−1Fn ext −Fn int+1 2(2−ch)˙ q0 a(5.15) El algoritmo de integraci´ on queda as´ ı completo a falta de determinar los par´ ametros de integraci´ on, como M,cyh. Para calcular cada componente diagonal de la matriz Mse utiliza la siguiente expre- si´ on Mii ≥h2 4 N X j=1 |Kij |(5.16) donde Kes la matriz global de rigidez, que en general es muy dif´ ıcil de obtener para sistemas at´ omicos. En su lugar, se podr´ ıa calcular la matriz cuasi-arm´ onica de rigidez de cada nodo como en [32], pero en este trabajo se ha optado por utilizar Mii =m, donde m es un par´ ametro que se estima con la ayuda de la matriz cuasi-arm´ onica de rigidez. Por otro lado, hes el paso de tiempo y normalmente toma el valor h= 1,0 y Nes el n´ umero de filas de la matriz K. El coeficiente de amortiguamiento se calcula como c≈2pλ0(5.17)
5. M´ etodos num´ ericos para la optimizaci´ on de problemas no lineales 26 Figura 5.2. Diagrama de bloques del m´ etodo del relajaci´ on din´ amica adaptativa.
5. M´ etodos num´ ericos para la optimizaci´ on de problemas no lineales 27 donde λ0es el menor autovalor de la matriz de Ky se calcula como λ0≈WT nSnWn WT nMnWn Sn ii =Ft int −Ft−1 int h˙ qt−1 2 i (5.18) donde Wes una matriz de peso elegida (en este caso la matriz identidad) y Sn ii es el componente ii de la matriz de rigidez direccional. En la Figura 5.2 se presenta un diagrama de bloques completo del m´ etodo implementado.
6. Potenciales termodin ´ amicos para el silicio A la hora de estudiar las propiedades de un material a nivel at´ omico, los resultados m´ as precisos se obtienen mediante c´ alculos ab initio, resolviendo la ecuaci´ on de Sch¨ ordinger, pero la complejidad es alta y el tama˜ no de la muestra que se puede analizar est´ a muy limitado. Alternativamente, se realizan simulaciones at´ omicas con m´ etodos como Monte Carlo o din´ amica molecular. Estos m´ etodos necesitan un modelo de fuerzas entre los ´ atomos que forman la muestra, los denominados potenciales interat´ omicos. Dada la importancia tecnol´ ogica del silicio, se han publicado y comparado [33, 34, 35] multitud de potenciales interat´ omicos para este material, con diferentes grados de sofisticaci´ on. A continuaci´ on se introducen dos de los m´ as usados: el potencial Stillinger- Weber [36] y el potencial del m´ etodo modificado del ´ atomo embebido [37], MEAM por sus siglas en ingl´ es. El primero tiene siete par´ ametros y la energ´ ıa total de una configuraci´ on at´ omica viene dada por t´ erminos de interacci´ on entre dos y tres ´ atomos. El segundo incluye enlaces direccionales y puede aplicarse a sistemas covalentes. Se presenta la formulaci´ on de ambos potenciales y los valores de los par´ ametros considerados en este estudio. Adem´ as, en los Ap´ endices A y B se recogen las expresiones de las fuerzas y del t´ ermino de expansi´ on por multipolos presentado anteriormente. 6.1. Potencial Stillinger-Weber La energ´ ıa potencial del sistema se aproxima como: E=1 2 N X i=1 X j,i vij(rij) + N X i=1 X j,iX k,i,j vjik(rij,rik) (6.1) donde el par´ ametro se elige de forma que f2toma como valor m´ ınimo −1, y σse elige para que el valor f221 6sea nulo. Adem´ as, f3posee simetr´ ıa traslacional y rotacional.
6. Potenciales termodin ´ amicos para el silicio 29 El termino de interacci´ on entre pares se escribe de la forma: vij(rij) = f2(rij/σ) (6.2) f2(r) = (A(Br−p−1)e(r−a)−1, r < a 0, r ≥a(6.3) donde los A,B,pyatoman valores positivos. Esta forma gen´ erica corta la funci´ on para r= asin discontinuidad en ninguna derivada, siendo esto una ventaja para las aplicaciones en din´ amica molecular. A continuaci´ on, en la Figura 6.1 se representa el valor de f2en funci´ on de rij/σ. Figura 6.1. Evaluaci´ on de la funci´ on f2(rij/σ). Por otro lado, el termino de interacci´ on de tres cuerpos toma la forma vjik(rij,rik) = f3(rij/σ,rik/σ) (6.4) f3(rij,rik) = h(rij,rik,θjik) + h(rji,rjk,θijk) + h(rki,rkj,θikj) (6.5) donde θjik es el ´ angulo entre rjyrksubtendido al v´ ertice i. La funci´ on hdepende de dos par´ ametros (λ,γ), y tiene la forma: h(rij,rik,θjik) = λeγ(rij−a)−1+γ(rik−a)−1cosjik +1 32 (6.6) El angulo tetrah´ edrico ideal θtes tal que cosθt=−1 3(6.7) de modo que la parte trigonom´ etrica de la expresi´ on discrimina a favor de pares de enlaces que salen del v´ ertice icon la geometr´ ıa deseada. Los valores utilizados para los par´ ametros que aparecen en las expresiones anteriores se presentan en la Tabla 6.1, donde se han modificado los valores originales de acuerdo con [38].
6. Potenciales termodin ´ amicos para el silicio 30 A 7.526106281 B 0.6022245584 p 4 a 1.8 λ22.4196 γ1.20 σ0.20951 2.1678 Tabla 6.1. Par´ ametros utilizados para el potencial Stillinger-Weber. 6.2. Potencial MEAM En el formalismo MEAM, la energ´ ıa total del sistema se aproxima como E=X i Fi(ρi) + 1 2X j,i φij(rij) (6.8) donde Fies la funci´ on embebida, ρies la densidad electr´ onica de fondo en el sitio i, y φij(rij) es el par de interacci´ on entre ´ atomos iyjseparados una distancia rij, respectivamente. De forma general, para calcular la energ´ ıa tienen que ser dadas las expresiones de las funciones Fiyφij(rij). La densidad electr´ onica de fondo se calcula teniendo en consideraci´ on la direccionalidad del enlace. La funci´ on embebida Fitiene una forma dada, pero no as´ ı el t´ ermino entre pares. En su lugar, se define una estructura de referencia donde los ´ atomos se encuentran en puntos exactos de la red y la energ´ ıa total por ´ atomo de dicha estructura se estima de la ecuaci´ on universal de estado a temperatura cero de Rose. Entonces, el t´ ermino entre pares se eval´ ua a partir de valores conocidos de la energ´ ıa total por ´ atomo y la energ´ ıa embebida, como una funci´ on de la distancia a sus vecinos. La energ´ ıa embebida tiene la forma Fi(ρi) = AiE0 iρilnρi(6.9) donde la energ´ ıa cohesiva E0 iy el par´ ametro Aidependen del tipo de ´ atomo. La densidad electr´ onica de fondo,ρiest´ a compuesta por una parte esf´ erica sim´ etrica, ρ(0) i, y otras partes angulares, ρ(1) i,ρ(2) i,ρ(3) i. Su expresi´ on se escribe como ρi=ρ(0) i ρ0Gi(Γi) (6.10) donde
6. Potenciales termodin ´ amicos para el silicio 31 Γi= 3 X k=1 t(k) i ρ(k) i ρ(0) i 2 G(Γ) = 2 1 + e−Γ (6.11) Las densidades electr´ onicas parciales vienen dadas por ρ(0) i=X j,i ρa(0) j(rij)Sij ρ(1) i= 3 X α=1 X j,i ρa(1) j(rij)rijα rij Sij 2 ρ(2) i= 3 X α=1 3 X β=1 X j,i ρa(2) j(rij)rijαrijβ r2 ij Sij 2 −1 3 X j,i ρa(2) j(rij)Sij 2 ρ(3) i= 3 X α=1 3 X β=1 3 X γ=1 X j,i ρa(3) j(rij)rijαrijβrijγ r3 ij Sij 2 −3 5 3 X α=1 X j,i ρa(3) j(rij)rijα rij Sij 2 (6.12) donde ρ(k) i(rij) = ρ0e"−β(k) i rij r0 ij −1!# (6.13) siendo ρ0un factor de escalado, β(h)ajustables y r0 ij la distancia al primer vecino en el equilibrio para la estructura de referencia. Como se coment´ o anteriormente, la funci´ on φij(rij) no se especifica directamente. En su lugar, viene dada por φij(rij) = φij(rij)Sij φij(rij) = 2 Zij0hEu i(rij)−Fi(ρ0(rij))i Eu i(rij) = −E0 ij(1 + a∗ ij +da∗3 ij )e−a∗ ij a∗ ij =α rij r0 ij −1 (6.14) donde Eu i(rij) es la funci´ on universal para una expansi´ on o contracci´ on uniforme en la estructura de referencia. En cuanto a la funci´ on de apantallamiento, est´ a dise˜ nada para que tome el valor Sij = 1 si los ´ atomos iyjno est´ an apantallados y dentro de un radio rc, mientras que
6. Potenciales termodin ´ amicos para el silicio 32 toma el valor Sij = 0 si est´ an completamente apantallados o fuera del radio de corte. La funci´ on es producto de una funci´ on de corte y una funci´ on de t´ erminos de tres cuerpos que considera el resto de ´ atomos del sistema Sij =Sijfcrc−rij ∆r Sij =Y k,i,j Sikj Sikj =fc Cikj −Cmin,ikj Cmax,ikj −Cmin,ikj ! Cikj = 1 + 2r2 ijr2 ik +r2 ijr2 jk −r4 ij r4 ij −(r2 ik −r2 jk)2 fc(x) = 1, x ≥1 h1−(1 −x)4i,0< x < 1 0, x ≤0 (6.15) El par´ ametro ∆rcontrola la distancia sobre la cual el radio de corte es suavizado desde 1 a 0 cerca de rij =rc. La versi´ on del potencial MEAM utilizada en este trabajo considera hasta segundos vecinos y los valores de los par´ ametros utilizados vienen dados en la Tabla 6.2. Ec4.63 β(2) 0.00 re2.35 β(3) 3.00 α4.90 t(1) 1.45 d 0.03 t(2) 7.61 A 0.53 t(3) -2.10 β(0) 3.00 Cmax 2.60 β(1) 7.5 Cmin 2.75 Tabla 6.2. Par´ ametros utilizados para el potencial MEAM. 6.3. Comparaci ´ on entre la aproximaci´ on de multipolos y la cuadratura de Gauss para el potencial Stillinger-Weber En este apartado se presenta una comparaci´ on entre las dos aproximaciones presentadas en el Cap´ ıtulo 3 para evaluar las medias de fase de campo medio. En concreto, se comparar´ an los resultados obtenidos a la hora de calcular distintas magnitudes con el potencial Stillinger-Weber, as´ ı como los tiempos de c´ alculo requeridos para ello. En las Figuras 6.2 y 6.3 se representa la energ´ ıa de formaci´ on, energ´ ıa por ´ atomo de la red perfecta, y el valor del par´ ametro de red de equilibrio para distintas temperaturas,
6. Potenciales termodin ´ amicos para el silicio 33 calculado minimizando el potencial Frespecto al valor del par´ ametro de red. Se comprueba que apenas hay diferencia entre ambas aproximaciones y, por tanto, ambas son v´ alidas para evaluar el potencial. Figura 6.2. Comparaci´ on de la energ´ ıa por ´ atomo de la red perfecta para distintas temperaturas, calculada con el potencial Stillinger-Weber y las dos aproximaciones presentadas: multipolos (MP) y cuadratura de Gauss (GAUSS). Figura 6.3. Comparaci´ on del par´ ametro de red de equilibrio para distintas temperaturas, calculado con el potencial Stillinger-Weber y las dos aproximaciones presentadas: multipolos (MP) y cuadratura de Gauss (GAUSS).
6. Potenciales termodin ´ amicos para el silicio 34 A continuaci´ on se presenta en la Tabla 6.3 la diferencia entre el valor de la energ´ ıa de formaci´ on calculada con ambas aproximaciones para distintas temperaturas y valores del par´ ametro de red. a/a0T=0K T=50K T=100K T=300K T=500K T=1000K 0.85 0.0 4 ·10−61,4·10−51,22 ·10−43,39 ·10−41,37 ·10−3 1.0 0.0 0.0 2,0·10−62,30 ·10−54,70 ·10−52,93 ·10−4 1.23 0.0 2,30 ·10−59,40 ·10−59,19 ·10−42,81 ·10−31,55 ·10−2 Tabla 6.3. Diferencia entre la energ´ ıa de formaci´ on, calculada por ambas aproximaciones para distintas temperaturas. De este modo, considerando que el valor proporcionado por la aproximaci´ on por multipolos es el real, el error m´ aximo que se da para los distintos valores del cociente a/a0y las distintas temperaturas ensayadas toma un valor de errmax = 0,696%, un valor que resulta admisible. En cuanto a las fuerzas, se ha utilizado un cubo con condiciones peri´ odicas de contorno con el ´ atomo que debe estar en el origen desplazado a la posici´ on [d, d, d] a distintas temperaturas y se ha calculado el m´ odulo de la fuerza del ´ atomo desplazado. En las Tablas 6.4, 6.5 y 6.6 se recogen las fuerzas, dadas en eV / ˚ A, calculadas utilizando cada una de las aproximaciones presentadas (MP y Gauss) de forma anal´ ıtica (A) y la diferencias entre ellas para distintos valores del desplazamiento d. Se observa como al aumentar el valor de dla diferencia entre las aproximaciones aumenta, as´ ı como al aumentar la temperatura para un valor de dfijo. T|FA |(MP) |FA |(Gauss) |FA(MP )−FA(Gauss)| 0 1.808792 1.808792 2,49 ·10−8 100 1.808902 1.808895 6,74 ·10−6 500 1.809342 1.809174 1,68 ·10−4 1000 1.809891 1.809241 6,5·10−4 Tabla 6.4. Fuerza de un ´ atomo calculada con las dos aproximaciones para distintas temperaturas y un valor de desplazamiento d= 0,01a. T|FA |(MP) |FA |(Gauss) |FA(MP )−FA(Gauss)| 0 0.167930 0.167930 1,16 ·10−13 100 0.167891 0.167890 6,00 ·10−7 500 0.167735 0.167720 1,54 ·10−5 1000 0.167540 0.167478 6,17 ·10−5 Tabla 6.5. Fuerza de un ´ atomo calculada con las dos aproximaciones para distintas temperaturas y un valor de desplazamiento d= 0,001a. En cuanto a los tiempos de c´ alculo, en las Tablas 6.7-6.9 se presentan los tiempos de c´ alculo de la energ´ ıa, la fuerza, donde FNies la fuerza calculada de forma num´ erica
6. Potenciales termodin ´ amicos para el silicio 41 La energ´ ıa potencial de la red Φse puede escribir utilizando un desarrollo de Taylor de la forma Φ=Φ0+X lj Φα l j!uα l j!+1 2X lj X l0j0 Φαβ l l0 j j0!uα l j!uβ l0 j0!+O(|u|3) (6.17) donde Φα l j!= ∂Φ ∂xαl j u=0 Φα l l0 j j0!= ∂2Φ ∂xαl j∂xβl0 j0 u=0 (6.18) La fuerza que act´ ua sobre el ´ atomo l jse puede calcular derivando Φcon respecto a u −F l j!=Φα l j!+X l0j0 Φαβ l l0 j j0!uβ l0 j0!+O(|u|2) (6.19) Teniendo en cuenta que en la configuraci´ on de equilibrio se da que Fl j= 0, la fuerza resulta F l j!=−X l0j0 Φαβ l l0 j j0!uβ l0 j0!(6.20) donde se han despreciado los t´ erminos de orden superior. Entonces, la matriz Φαβl l0 j j0 representa la fuerza cambiada de signo que aparece en el ´ atomo l jen la direcci´ on α cuando se aplica un desplazamiento unitario en el ´ atomo l0 j0en la direcci´ on β. Imponiendo las distintas condiciones de simetr´ ıa de la red es posible establecer relaciones entre las matrices de los ´ atomos que se encuentran a la misma distancia del ´ atomo de referencia conformando los distintos niveles de vecinos. En el caso de este trabajo se va a considerar hasta segundos vecinos. Las matrices de constantes de fuerza de primeros vecinos para la estructura c´ ubica de diamante tienen, de acuerdo con [39], la siguente forma α β β β α β β β α α−β−β −β α β −β β α α−β β −β α −β β−β α α β −β β α −β −β−β α Por otro lado, las matrices correspondientes a segundos vecinos tienen la forma
6. Potenciales termodin ´ amicos para el silicio 42 µ ν δ ν µ δ γγλ µ−ν δ −ν µ −δ γ−γ λ µ−ν−δ −ν µ δ −γ γ λ µ ν −δ ν µ −δ −γ−γ λ µ δ ν γ λ γ ν δ µ µ δ −ν γ λ −γ −ν−δ µ µ−δ−ν −γ λ γ −ν δ µ µ−δ ν −γ λ −γ ν−δ µ λ γ γ δ µ ν δ ν µ λ γ −γ δ µ −ν −δ−ν µ λ−γ γ −δ µ −ν δ−ν µ λ−γ−γ −δ µ ν −δ ν µ En este trabajo se han calculado las 7 constantes (α,β,µ,ν,δ,γyλ) utilizando los dos potenciales presentados. En la Tabla 6.14 se comparan los resultados obtenidos con los valores calculados mediante primeros principios [40]. Constante SW MEAM Ab initio α-7.8115 -5.4678 -5.47 β-4.9131 -4.1257 -3.88 µ-0.1207 -0.368 -0.35 ν-0.1207 -0.3128 -0.29 δ0.2415 0.3776 0.68 γ-0.2415 -0.3776 -0.68 λ0.4831 0.0103 0.17 Tabla 6.14. Comparaci´ on de los valores de las constantes de fuerza para primeros y segundos vecinos del silicio calculadas con distintos m´ etodos. Los valores est´ an dados en 104dyn/cm. Se puede comprobar como el potencial MEAM es el que ofrece unos resultados m´ as cercanos a los obtenidos mediante primeros principios. Adem´ as, se han calculado las constantes el´ asticas del silicio, utilizando las constantes de fuerza calculadas, para comparar tambi´ en la validez de ambos potenciales. En concreto, las constantes el´ asticas se pueden calcular a partir de las constantes de fuerzas mediante las expresiones de la Ecuaci´ on 6.21 c11 =−(α+ 8µ)/a c12 = (α−2β+ 4µ−8ν+ 4λ)/a c44 = (−α−4µ−4λ+β2/α)/a (6.21) donde aes el par´ ametro de red y toma el valor ya mencionado. A continuaci´ on, se presenta en la Tabla 6.15 una comparaci´ on entre los valores de las constantes el´ asticas obtenidas con los potenciales presentados, con el m´ etodo de primeros principios usado en la comparaci´ on anterior y con los valores experimentales [41]
6. Potenciales termodin ´ amicos para el silicio 43 Constante SW MEAM Ab initio Experimental c11 161.6 154.9 152.3 165.7 c12 81.6 71.0 71.6 63.9 c44 60.3 69.7 63.3 79.6 Tabla 6.15. Comparaci´ on de los valores de las constantes el´ asticas del silicio calculadas con distintos m´ etodos. Los valores est´ an dados en GP a. Se observa como de nuevo el potencial MEAM presenta la mejor aproximaci´ on a los valores experimentales. 6.5.2. Matriz din ´ amica Utilizando la matriz de constantes de fuerza definida en la secci´ on la ecuaci´ on de movimiento del ´ atomo l jse puede escribir como mj¨ uα l j!+X l0j0 Φαβ l−l0 j j0!uβ l0 j0!= 0 (6.22) donde mjes la masa del ´ atomo jen la celda. Considerando la variable reducida wα(j) = m1/2 juα(j) (6.23) la Ecuaci´ on 6.22 resulta ¨ wα(j) + X l0j0 Dαβ l−l0 j j0!wβ l0 j0!= 0 (6.24) donde Dαβ l−l0 j j0!=1 √mjmj0 Φαβ l−l0 j j0!(6.25) es la matriz din´ amica. Ensayando una soluci´ on del tipo wα l j!=Wα(j)eiωte−ik·r(l j)(6.26) se tiene
6. Potenciales termodin ´ amicos para el silicio 44 ω2Wα l j!=X l0j0 Dαβ l−l0 j j0!Wβ(j0)e−ik·hr(l j)−r(l0 j0)i(6.27) resultando un problema de autovalores y autovectores. En el silicio la red est´ a compuesta por dos tipos de ´ atomos, uno por cada una de las dos redes FCC que se entrecruzan, de modo que la matriz din´ amica tiene la expresi´ on D(k) = "DAA(k)DAB(k) DBA(k)DAA(k)# Resolviendo el problema de autovalores y autovectores planteado en la Ecuaci´ on 6.27 se obtienen las frecuencias de vibraci´ on para distintos valores del vector de onda kpara cada uno de los 6 modos de vibraci´ on, 3 por cada tipo de ´ atomo. En el caso del silicio, considerando que el ´ atomo de referencia es de tipo A, todos los primeros vecinos ser´ an del tipo By todos los segundos vecinos ser´ an del tipo A. De este modo, las submatrices que forman la matriz din´ amica ser´ an de la forma DAA(k) = Φ0+X i∈Ω1 Φie−ik·ri DAB(k) = DBA(k) = X i∈Ω2 Φie−ik·ri (6.28) donde Ω1yΩ2son los conjuntos de primeros y segundos vecinos respectivamente. 6.5.3. Dispersi ´ on de fonones La relaci´ on ω(k) se conoce como relaci´ on de dispersi´ on y es muy ´ util a la hora de estudiar procesos de transporte puesto que contiene informaci´ on de la forma en la que vibra la red. Debido a la periodicidad de la red, esta relaci´ on solo es necesaria calcularla dentro de una peque˜ na zona conocida como Primera Zona de Brillouin, de modo que los valores del vector de onda que se estudian son los contenidos en esta zona. Adem´ as, se estudian normalmente direcciones de alta simetr´ ıa, que vienen definidas generalmente por una serie de puntos caracter´ ısticos. La Primera Zona de Brillouin para el caso del silicio se muestra en la Figura 6.8. Las posiciones de los puntos de alta simetr´ ıa se˜ nalados en la figura se presentan a continuaci´ on Γ=2π a 0 0 0 X=2π a 1 0 0 L=2π a 1 2 1 2 1 2
6. Potenciales termodin ´ amicos para el silicio 45 W=2π a 1 1 2 0 K=2π a 3 4 3 4 0 U=2π a 1 1 4 1 4 Figura 6.8. Primera zona de Brillouin en la estructura c´ ubica del diamante. En esta estructura, las lineas de simetr´ ıa a lo largo de las cuales se estudia la disper- si´ on son estas tres ∆:k=2π a(ξ,0,0) 0 < ξ < 1 Λ:k=2π a(ξ,ξ,ξ) 0 < ξ < 1 2 Σ:k=2π a(ξ,ξ,0) 0 < ξ < 1 (6.29) De esta forma se han obtenido las curvas de fonones para ambos potenciales. Adem´ as, para el caso del potencial Stillinger-Weber se han considerado dos modelos: uno hasta primeros vecinos y otro hasta segundos. En las Figuras 6.9, 6.10 y 6.11 se comparan los resultados obtenidos con resultados experimentales [42, 43]. Se observa como ambos potenciales captan bien la tendencia y las magnitudes observadas en los experimentos, aunque presentan mayor r´ ıgidez.
6. Potenciales termodin ´ amicos para el silicio 46 Figura 6.9. Curva de dispersi´ on de fonones obtenida con un modelo hasta primeros vecinos del potencial Stillinger-Weber (l´ ınea continua) frente a resultados experimentales (c´ ırculos [43], tri´ angulos [42]). Figura 6.10. Curva de dispersi´ on de fonones obtenida con un modelo hasta segundos vecinos del potencial Stillinger-Weber (l´ ınea continua) frente a resultados experimentales (c´ ırculos [43], tri´ angulos [42]).
6. Potenciales termodin ´ amicos para el silicio 47 Figura 6.11. Curva de dispersi´ on de fonones obtenida con un modelo hasta segundos vecinos del potencial MEAM (l´ ınea continua) frente a resultados experimentales (c´ ırculos [43], tri´ angulos [42]). Figura 6.12. Curva de dispersi´ on de fonones obtenida con un modelo hasta primeros vecinos del potencial Stillinger-Weber (azul) frente a resultados obtenidos en [33] (amarillo).
6. Potenciales termodin ´ amicos para el silicio 48 Figura 6.13. Curvas de dispersi´ on de fonones obtenida con un modelo hasta segundos vecinos del potencial Stillinger-Weber (azul) y del potencial MEAM (rojo) frente a resultados obtenidos en [33] (amarillo). Por ´ ultimo se comparan los resultados obtenidos en este estudio con otros resultados computacionales obtenidos con din´ amica de red en [44], utilizando el potencial EDIP (Environment Dependent Interatomic Potential)[33]. En la Figura 6.12 se comparan los resultados obtenidos con el modelo hasta primeros vecinos del potencial Stillinger-Weber con los de la fuente citada, mientras que en la Figura 6.13 se comparan los resultados del modelo hasta segundos vecinos de los potenciales utilizados en este estudio con los de la fuente. Se comprueba as´ ı que los resultados obtenidos con la misma metodolog´ ıa aqu´ ı empleada y distinto potencial presentan resultados similares a los presentados anteriormente.
7. Conclusiones y trabajos futuros En este trabajo se ha aplicado una teor´ ıa termodin´ amica estad´ ıstica de no equilibrio para el estudio procesos termomec´ anicos en materiales a escala nanosc´ opica. Se ha presentado un modelo simple para calcular la conductividad t´ ermica en nanohilos de silicio, se han ajustado los coeficientes emp´ ıricos de dicho modelo y se han comparado los resultados del mismo con otros resultados experimentales, quedando demostrada la capacidad del mismo para reproducir el comportamiento de este nanomaterial. Tambi´ en se han presentado dos potenciales interat´ omicos de amplia utilizaci´ on para el estudio termomec´ anico del silicio y dos m´ etodos num´ ericos para la resoluci´ on de sistemas de ecuaciones no lineales que se utilizan para resolver los sistemas que aparecen en la teor´ ıa comentada. Para la validaci´ on de la teor´ ıa y los potenciales presentados se ha obtenido distintas magnitudes f ˜ Asicas como la energ´ ıa de formaci´ on, las constantes de fuerza, las constantes el´ asticas o las curvas de dispersi´ on de fonones y se han comparado ambos potenciales entre s´ ı y con otros resultados tanto experimentales como num´ ericos, poni´ endose de manifiesto la idoneidad tanto de la teor´ ıa como de los potenciales para su utilizaci´ on en procesos termomec´ anicos en silicio. Los posibles trabajos futuros se centran principalmente en las siguientes l´ ıneas: Construcci´ on de un modelo m´ as complejo para el estudio de conductividad t´ ermica en nanohilos de silicio que involucre la utilizaci´ on de los potenciales presentados. Estudio de la influencia de defectos en las propiedades termomec´ anicas de nanohilos de silicio. Estudio del comportamiento de otras estructuras como las nanopart´ ıculas de silicio. Aplicaci´ on de la teor´ ıa presentada a otros materiales de inter´ es tecnol´ ogico y actualidad como el grafeno.
A. Ap´ endice A A.1. Derivadas del potencial Stillinger-Weber Conociendo la expresi´ on de la energ´ ıa para el potencial Stillinger-Weber E=1 2 N X i=1 X j,i vij(rij) + N X i=1 X j,iX k,i,j vjik(rij,rik) (A.1) , la fuerza se puede escribir entonces como fi=∂E ∂ri =−1 2 X j,i ∂vij(rij) ∂ri +X j,i ∂vij(rji) ∂ri− +X j,iX k,i,j ∂vjik(rij,rik) ∂ri + 2 ·X j,iX k,i,j ∂vjik(rji,rjk) ∂ri (A.2) A.1.1. Derivadas del t´ ermino vij(rij) ∂vij(rij) ∂ri =∂vij(rij) ∂rij ∂rij ∂rij ∂rij ∂ri =−∂vij(rij) ∂rij rij rij (A.3) ∂vij(rij) ∂rij =−A σerij σ−a−1"pB σrij σ−p−1 + +Brij σ−p −1 1 rij σ−a 2 (A.4)
B. Ap´ endice B 57 B.1.3. Derivadas del t´ ermino ρ(k) ij Y1iσ =X j,i ρa(1) j rijσ rij Sij (B.25) Y2iσβ =X j,i ρa(2) j rijσ rijβ r2 ij Sij (B.26) Y3iσβγ =X j,i ρa(3) j rijσ rijβrijγ r3 ij Sij (B.27) W2i=X j,i ρa(2) jSij (B.28) W3iσ =X j,i ρa(3) j rijσ rij Sij (B.29) ρ(1) i2 = 3 X σ=1 Y2 1iσ (B.30) ρ(2) i2 = 3 X σ=1 3 X β=1 Y2 2iσβ −1 3W2 2i(B.31) ρ(3) i2 = 3 X σ=1 3 X β=1 3 X γ=1 Y2 3iσβγ −3 5 3 X σ=1 W2 3iσ (B.32) ∂ρ(k) i2 ∂rlm =∂ρ(k)r O im δil +∂ρ(k)r O il δim (B.33) ∂ρ(k) i2 ∂Slm =∂ρ(k)S O im δil +∂ρ(k)S O il δim (B.34) ∂ρ(k) i2 ∂rlm =∂ρ(k) imαδil −∂ρ(k) ilαδim (B.35) ∂ρ(1) r O ij =2 rij Sij ρa(1)0 j(rij)−1 rij ρa(1) j(rij)! X σ Y1iσ rijσ (B.36)
B. Ap´ endice B 58 ∂ρ(1) S O ij =2 rij ρa(1) j(rij)X σ Y1iσ rijσ (B.37) ∂ρ(1) ijα =2 rij SijY1iαρa(1) j(rij) (B.38) ∂ρ(2) r O ij =2 r2 ij Sij ρa(2)0 j(rij)−2 rij ρa(2) j(rij)! X σX β Y2iσβrijσ rijβ −2 3W3iσ ρa(2)0 j(rij)Sij (B.39) ∂ρ(2) S O ij =2 r2 ij ρa(2) j(rij) X σX β Y2iσβrijσ rijβ −2 3W3iσ ρa(2) j(rij)Sij (B.40) ∂ρ(2) ijα =4 r2 ij ρa(2) j(rij)Sij X σ Y2iσβrijσ (B.41) ∂ρ(3) r O ij =2 r3 ij Sij ρa(3)0 j(rij)−3 rij ρa(3) j(rij)! X σX βX γ Y3iσβγrijσ rijβrijγ − −6 5rij Sij ρa(3)0 j(rij)−1 rij ρa(3) j(rij)! X σ W3iσ rijσ (B.42) ∂ρ(3) S O ij =2 r3 ij ρa(3) j(rij)X σX βX γ Y3iσβγrijσ rijβrijγ− −6 5rij ρa(3) j(rij)X σ W3iσ (B.43) ∂ρ(3) ijα =6 r3 ij ρa(3) j(rij)Sij X σX β Y3iασβrijσ rijβ −6 5rij Sijρa(3) j(rij)W3iα (B.44) B.1.4. Derivadas del t´ ermino t(k) i ∂t(k) i ∂rlm =∂t(k)r O im δil +∂t(k)r O il δim (B.45) ∂t(k) i ∂Sij =∂t(k)S O im δil +∂t(k)S O il δim (B.46)
B. Ap´ endice B 59 ∂t(k)r O ij =1 ρ(0) it(k) 0,j −t(k) iρa(0)0 j(rij)Sij (B.47) ∂t(k)S O ij =1 ρ(0) it(k) 0,j −t(k) iρa(0) j(rij) (B.48) B.1.5. Derivadas del t´ ermino Γ(k) i ∂Γ(k) i ∂rlm =∂Γr O im δil +∂Γr O il δim (B.49) ∂Γ(k) i ∂Sij =∂ΓS O im δil +∂ΓS O il δim (B.50) ∂Γ(k) i ∂rlmα =∂Γimαδil +∂Γilαδim (B.51) ∂Γr O ij =−2 ρ(0) i Γi∂ρ(0) r O ij +ρa(0)0 j(rij)Sij ρ(0) i3 3 X k=1 t(k) 0,j −t(k) iρ(k) i2 + +1 (ρ(0) ij )2 3 X k=1 t(k) i∂ρ(k)S O ij (B.52) Γijα =1 (ρ(0) ij )2 3 X k=1 t(k) i∂ρ(k) ij (B.53) B.1.6. Derivadas del t´ ermino ρi ∂ρ(k) i ∂rlm =∂ρ r O imδil +∂ρ r O il δim (B.54) ∂ρ(k) i ∂Sij =∂ρ S O imδil +∂ρ S O il δim (B.55) ∂ρ(k) i ∂rlmα =∂ρimαδil +∂ρilαδim (B.56)
B. Ap´ endice B 60 ∂ρ r O ij =1 ρ0 i(Gi(Γi)−2G0 i(Γi)Γi)∂ρ(0) r O ij + +1 ρ(0) i G0 i(Γi) 3 X k=1 ∂t(k)r O ij (ρ(k) i)2+t(k) i∂ρ(k)r O ij − −ρiρi0 Zi0 G0 i(Γref i) 3 X k=1 s(k) i∂t(k)r O ij (B.57) ∂ρ S O ij =1 ρ0 i(Gi(Γi)−2G0 i(Γi)Γi)∂ρ(0) S O ij + +1 ρ(0) i G0 i(Γi) 3 X k=1 ∂t(k)S O ij (ρ(k) i)2+t(k) i∂ρ(k)S O ij − −ρiρi0 Zi0 G0 i(Γref i) 3 X k=1 s(k) i∂t(k)S O ij (B.58) ∂ρijα =1 ρ0 iρ(0) i G0 i(Γi) 3 X k=1 t(k) i∂ρ(k) ijα (B.59) B.2. Expresi ´ on del t´ ermino ∆jVi(q)para el MEAM En este caso, debido a la complejidad del potencial, el t´ ermino se calcular´ a de forma num´ erica, utilizando diferencias centrales de segundo orden f00(x0)≈f(x0+h)−2f(x0) + f(x0−h) h2(B.60) De este modo, el t´ ermino se aproxima como ∆lVi({q}) = 3 X α=1 ∂2 qlα Vi({q}) = 3 X α=1 ∂2Vi ∂q2 lα ≈ 3 X α=1 Vi(qi,qj,...,qlα +h,qN)−Vi({q}) + Vi(qi,qj,...,qlα −h,qN) h2 (B.61)
´ Indice de figuras 2.1. Estructura cristalina c´ ubica de diamante, presente en el silicio. . . . . . . . 5 2.2. Conductividad t´ ermica del silicio macizo frente a la temperatura. . . . . . . 5 2.3. Conductividad t´ ermica de distintos nanohilos de silicio frente a la temperatura......................................... 6 5.1. Diagrama de bloques del m´ etodo del gradiente conjugado no lineal. . . . . 23 5.2. Diagrama de bloques del m´ etodo del relajaci´ on din´ amica adaptativa. . . . 26 6.1. Evaluaci´ on de la funci´ on f2(rij/σ)......................... 29 6.2. Comparaci´ on de la energ´ ıa por ´ atomo de la red perfecta para distintas temperaturas, calculada con el potencial Stillinger-Weber y las dos aproximaciones presentadas: multipolos (MP) y cuadratura de Gauss (GAUSS). . 33 6.3. Comparaci´ on del par´ ametro de red de equilibrio para distintas temperaturas, calculado con el potencial Stillinger-Weber y las dos aproximaciones presentadas: multipolos (MP) y cuadratura de Gauss (GAUSS). . . . . . . . 33 6.4. Energ´ ıa por ´ atomo de la red perfecta calculada con los potenciales SW y MEAM para distintas temperaturas usando la aproximaci´ on multipolos. . . 36 6.5. Valor del funcional Fpor ´ atomo de la red perfecta calculada con los potenciales SW y MEAM para distintas temperaturas usando la aproximaci´ on multipolos...................................... 36 6.6. Energ´ ıa por ´ atomo de la red para distintos valores del par´ ametro de red y distintas temperaturas, calculada con los potenciales SW y MEAM utilizando la aproximaci´ onmultipolos. ....................... 38 i
Bibliograf´ ıa ii 6.7. Valor del funcional Fpor ´ atomo de la red para distintos valores del par´ ametro de red y distintas temperaturas, calculada con los potenciales SW y MEAM utilizando la aproximaci´ on multipolos. . . . . . . . . . . . . . . . . 39 6.8. Primera zona de Brillouin en la estructura c´ ubica del diamante. . . . . . . . 45 6.9. Curva de dispersi´ on de fonones obtenida con un modelo hasta primeros vecinos del potencial Stillinger-Weber (l´ ınea continua) frente a resultados experimentales (c´ ırculos [43], tri´ angulos [42]). . . . . . . . . . . . . . . . . . 46 6.10.Curva de dispersi´ on de fonones obtenida con un modelo hasta segundos vecinos del potencial Stillinger-Weber (l´ ınea continua) frente a resultados experimentales (c´ ırculos [43], tri´ angulos [42]). . . . . . . . . . . . . . . . . . 46 6.11.Curva de dispersi´ on de fonones obtenida con un modelo hasta segundos vecinos del potencial MEAM (l´ ınea continua) frente a resultados experimentales (c´ ırculos [43], tri´ angulos[42]). .................... 47 6.12.Curva de dispersi´ on de fonones obtenida con un modelo hasta primeros vecinos del potencial Stillinger-Weber (azul) frente a resultados obtenidos en[33](amarillo).................................. 47 6.13.Curvas de dispersi´ on de fonones obtenida con un modelo hasta segundos vecinos del potencial Stillinger-Weber (azul) y del potencial MEAM (rojo) frente a resultados obtenidos en [33] (amarillo). . . . . . . . . . . . . . . . . 48
Bibliograf´ ıa [1] David Chandler. Introduction to modern statistical mechanics. Introduction to Modern Statistical Mechanics, by David Chandler, pp. 288. Foreword by David Chandler. Oxford University Press, Sep 1987. ISBN-10: 0195042778. ISBN-13: 9780195042771, 1, 1987. [2] John Hammersley. Monte carlo methods. Springer Science & Business Media, 2013. [3] Wm G Hoover. Molecular dynamics. In Molecular Dynamics, volume 258, 1986. [4] G. Venturini, K. Wang, I. Romero, M.P. P Ariza, and M. Ortiz. Atomistic long-term simulation of heat and mass transport. Journal of the Mechanics and Physics of Solids, 73:242–268, 2014. [5] L Weber and E Gmelin. Transport properties of silicon. Applied Physics A, 53(2):136– 140, 1991. [6] C. J. Glassbrenner and Glen A. Slack. Thermal conductivity of silicon and germanium from 3??K to the melting point. Physical Review, 1964. [7] Glen A. Slack. Thermal conductivity of pure and impure silicon, silicon carbide, and diamond. Journal of Applied Physics, 1964. [8] Allon I Hochbaum, Renkun Chen, Raul Diaz Delgado, Wenjie Liang, Erik C Garnett, Mark Najarian, Arun Majumdar, and Peidong Yang. Enhanced thermoelectric performance of rough silicon nanowires. Nature, 451(January):163–167, 2008. [9] Peidong Yang Peidong Yang. Semiconductor nanowires for energy conversion. Nanoelectronics Conference (INEC), 2010 3rd International, pages 527–546, 2010. [10] Akram I Boukai, Yuri Bunimovich, Jamil Tahir-Kheli, Jen-Kan Yu, William A Goddard, and James R Heath. Silicon nanowires as efficient thermoelectric materials. Nature, 451(7175):168–71, 2008. [11] C Chiappini, E De Rosa, J O Martinez, X Liu, J Steele, M M Stevens, and E Tasciotti. Biodegradable silicon nanoneedles delivering nucleic acids intracellularly induce localized in vivo neovascularization. Nature materials, 14(March):6–13, 2015. iii
Bibliograf´ ıa iv [12] Elizabeth Buitrago, Montserrat Fern´ andez Bola˜ nos Badia, Yordan M. Georgiev, Ran Yu, Olan Lotty, Justin D. Holmes, Adrian M. Nightingale, Ho¨ el M. Guerin, and Adrian M. Ionescu. Electrical characterization of high performance, liquid gated vertically stacked SiNW-based 3D FET biosensors. Sensors and Actuators, B: Chemical, 199:291–300, 2014. [13] Guo-Jun Zhang, Li Zhang, Min Joon Huang, Zhan Hong Henry Luo, Guang Kai Ignatius Tay, Eu-Jin Andy Lim, Tae Goo Kang, and Yu Chen. Silicon nanowire biosensor for highly sensitive and rapid detection of Dengue virus. Sensors and Actuators B: Chemical, 146(1):138–144, 2010. [14] C K Chan, H Peng, G Liu, K McIlwrath, X F Zhang, R A Huggins, and Y Cui. High-performance lithium battery anodes using silicon nanowires. Nat Nanotechnol, 3(1):31–35, 2008. [15] John M Ziman. Electrons and phonons: the theory of transport phenomena in solids. Oxford University Press, 1960. [16] Deyu Li, Yiying Wu, Philip Kim, Li Shi, Peidong Yang, and Arun Majumdar. Thermal conductivity of individual silicon nanowires. Applied Physics Letters, 83(14):2934–2936, 2003. [17] L. Onsager. Irreversible processes. Phys. Rev., 37:237–241, 1931. [18] Lars Onsager. Reciprocal relations in irreversible processes. I. Physical Review, 37(4):405–426, 1931. [19] Julia M Yeomans. Statistical mechanics of phase transitions. Clarendon Press, 1992. [20] Arthur H Stroud. Approximate calculation of multiple integrals. 1971. [21] Yong-Hee Park, Jungwon Kim, Hyoungjoon Kim, Ilsoo Kim, Ki-Young Lee, Dongjea Seo, Heon-Jin Choi, and Woochul Kim. Thermal conductivity of vls-grown rough si nanowires with various surface roughnesses and diameters. Applied Physics A, 104(1):7–14, 2011. [22] Etienne Blandre, Laurent Chaput, Samy Merabia, David Lacroix, and Konstantinos Termentzidis. Modeling the reduction of thermal conductivity in core/shell and diameter-modulated silicon nanowires. Physical Review B, 91(11):115404, 2015. [23] Li-Feng Cui, Riccardo Ruffo, Candace K Chan, Hailin Peng, and Yi Cui. Crystallineamorphous coreshell silicon nanowires for high capacity and high current battery electrodes. Nano Letters, 9(1):491–495, 2008. [24] Pierre Martin, Zlatan Aksamija, Eric Pop, and Umberto Ravaioli. Impact of phononsurface roughness scattering on thermal conductivity of thin si nanowires. Physical review letters, 102(12):125503, 2009. [25] Frederic Sansoz. Surface faceting dependence of thermal transport in silicon nanowires. Nano letters, 11(12):5378–5382, 2011. [26] FM Ross, J Tersoff, and MC Reuter. Sawtooth faceting in silicon nanowires. Physical review letters, 95(14):146104, 2005.
Bibliograf´ ıa v [27] David G Cahill, Henry E Fischer, Tom Klitsner, ET Swartz, and RO Pohl. Thermal conductivity of thin films: measurements and understanding. Journal of Vacuum Science & Technology A, 7(3):1259–1266, 1989. [28] Carl T Kelley. Solving nonlinear equations with Newton’s method, volume 1. Siam, 2003. [29] Jonathan Richard Shewchuk. An Introduction to the Conjugate Gradient Method Without the Agonizing Pain. 1994. [30] Youcef Saad and Martin H Schultz. Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems. SIAM Journal on scientific and statistical computing, 7(3):856–869, 1986. [31] Charles G Broyden. A class of methods for solving nonlinear simultaneous equations. Mathematics of computation, 19(92):577–593, 1965. [32] Yashashree Kulkarni. Coarse-graining of atomistic description at finite temperature. PhD thesis, California Institute of Technology, 2007. [33] Martin Z Bazant, Efthimios Kaxiras, and JF Justo. Environment-dependent interatomic potential for bulk silicon. Physical Review B, 56(14):8542, 1997. [34] Stephen J Cook and Paulette Clancy. Comparison of semi-empirical potential functions for silicon and germanium. Physical Review B, 47(13):7686, 1993. [35] Carolina Abs da Cruz, Konstantinos Termentzidis, Patrice Chantrenne, and Xavier Kleber. Molecular dynamics simulations for the prediction of thermal conductivity of bulk silicon and silicon nanowires: Influence of interatomic potentials and boundary conditions. Journal of Applied Physics, 110(3):034309, 2011. [36] Frank H. Stillinger and Thomas A. Weber. Computer simulation of local order in condensed phases of silicon. Physical Review B, 1985. [37] Byeong-Joo Lee. A modified embedded atom method interatomic potential for silicon. Calphad, 31(1):95–104, 2007. [38] BJ Thijsse. Relationship between the modified embedded-atom method and stillinger-weber potentials in calculating the structure of silicon. Physical Review B, 65(19):195207, 2002. [39] Frank Herman. Lattice vibrational spectrum of germanium. Journal of Physics and Chemistry of Solids, 8:405–418, 1959. [40] Siqing Wei and MY Chou. Ab initio calculation of force constants and full phonon dispersions. Physical review letters, 69(19):2799, 1992. [41] Walter A Harrison. Electronic structure and the properties of solids: the physics of the chemical bond. Courier Corporation, 2012. [42] G Nilsson and G Nelin. Study of the homology between silicon and germanium by thermal-neutron spectrometry. Physical Review B, 6(10):3777, 1972. [43] G Dolling. Lattice vibrations in crystals with the diamond structure. In Inelastic Scattering of Neutrons in Solids and Liquids. V. II. Proceedings of the Symposium on Inelastic Scattering of Neurons in Solids and Liquids, 1963.
Bibliograf´ ıa vi [44] Asegun S Henry and Gang Chen. Spectral phonon transport properties of silicon based on molecular dynamics simulations and lattice dynamics. Journal of Computational and Theoretical Nanoscience, 5(2):141–152, 2008.