Repositorio Institucional de Documentos

UNIVERSIDAD DE ZARAGOZA FACULTAD DE CIENCIAS DEPARTAMENTO DE F´ ISICA DE LA MATERIA CONDENSADA TRABAJO DE FINAL DE GRADO Propiedades optoelectr´onicas del grafeno ´ I˜nigo Arricibita Yoldi Tutor: Luis Mart´ın Moreno ´ Indice 1. Introducci´on 2 1.1. Objetivos .............................................. 2 1.2. Fundamentos ............................................ 3 1.2.1. Propagaci´on de ondas electromagn´eticas en el vac´ıo . . . . . . . . . . . . . . . . . . 3 1.2.2. C´alculo de coeficientes de transmisi´on y reflexi´on, confinamiento de modos p.... 5 2. Modelo 8 2.1. Planteamiento: amplitud como ecuaci´on integral . . . . . . . . . . . . . . . . . . . . . . . . 8 2.2. Ejemplo de aplicaci´on: defecto gaussiano . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 3. Perfil de impurezas de Ngaussianas 13 3.1. Aplicaci´ondelmodelo ....................................... 13 3.2. Predicci´onenFOBA........................................ 15 3.3. Resultados de las simulaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3.1. Discretizaci´on de la ecuaci´on integral . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 3.3.2. Comparaci´on FOBA anal´ıtica y de simulaciones . . . . . . . . . . . . . . . . . . . . . 19 3.3.3. Comparaci´on FOBA y simulaciones . . . . . . . . . . . . . . . . . . . . . . . . . . . 20 3.3.4. Comportamiento con a,anchuradeldefecto....................... 21 3.3.5. Comportamiento con δ,alturadeldefecto........................ 22 3.3.6. Comportamiento con N, n´umero de gaussianas . . . . . . . . . . . . . . . . . . . . . 23 3.3.7. Comportamiento con γ, distancia entre gaussianas . . . . . . . . . . . . . . . . . . . 25 4. Conclusiones 26 5. Ap´endice 27 5.1. Propiedades de la base de modos TM y TE . . . . . . . . . . . . . . . . . . . . . . . . . . . 27 5.2. C´alculo de la transformada de fourier en el defecto Gaussiano . . . . . . . . . . . . . . . . . 29 5.3. C´alculo de coeficiente de reflexi´on R............................... 30 5.4. C´alculo de integrales con la funci´on de Green . . . . . . . . . . . . . . . . . . . . . . . . . . 31 5.5. C´alculo de la integral de una gaussiana . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32 5.6. C´alculo de la suma geom´etrica . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 5.7. C´alculo del m´aximo de reflexi´on en FOBA . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 5.8. C´odigoMatLab........................................... 34 1. Introducci´on 1.1. Objetivos El grafeno es un material formado por una capa monoat´omica de carbono formando una extructura hexagonal peri´odica. Su origen se remonta al grafito (cuya estructura se resolvi´o en 1916 [1]), que consta, esencialmente, de varias capas de grafeno; aunque ´este t´ermino no se toma hasta 1987 [3]. En 1949, Philip Russel Wallace calcul´o la estructura de bandas de este material [2]. Seg´un estos c´alculos, se predec´ıa que el grafeno era una estructura inestable y es por eso que se tard´o en comenzar a tratar de obtener una sola l´amina de grafeno. A nivel experimental, se comenz´o trabajando con l´aminas de grafito muy finas: las primeras im´agenes (tomadas por microscopio electr´onico de transmisi´on) de grafito de pocas capas datan de 1948 [4]. Posteriormente, se lleg´o a detectar grafito del grosor de un ´atomo [5], lo que condujo al crecimiento epitaxial de grafeno en otros materiales [6]. Estos resultados de 1997, no obstante, no daban como resultado grafeno, pues para considerarlo como tal deb´ıa estar en vac´ıo y no crecido sobre otro material (ya que entonces se produce hibridaci´on entre los orbitales del grafeno y los del material sobre el que se crece). No fue hasta 2004 cuando se consigui´o aislar una capa de grafeno aislada. Andre Geim y Kostya Novoselov, de la universidad de Manchester, lograron aislar capas de grafeno a trav´es de grafito mediante la t´ecnica de ((cinta adhesiva Scotch)) [7]. En 2010 se les otorg´o el premio Nobel por este trabajo. Durante los ´ultimos a˜nos se ha estudiado el grafeno debido a sus m´ultiples propiedades, como por ejemplo su flexibilidad y elasticidad [8] o sus altas conductividades el´ectrica [9] y t´ermica [10]. En este trabajo nos centraremos en su capacidad como conductor de ondas electromagn´eticas: desarrollaremos un modelo de propagaci´on de plasmones de superficie en grafeno (graphene surface plasmons, GSP), es decir, c´omo ciertos modos de ondas electromagn´eticas (transversales magn´eticos) pueden confinarse en la superficie del grafeno y c´omo, a trav´es de impurezas en la conductividad (bien inducidas por un potencial o bien propias del material) se puede modificar la propagaci´on de tales plasmones. Concretamente, el estudio de este trabajo se focaliza en unas impurezas distribuidas espacialmente en forma de gaussianas para la conductividad el´ectrica. Al final del trabajo, discutiremos las posibles aplicaciones de este sistema para casos experimentales. 2 1.2. Fundamentos 1.2.1. Propagaci´on de ondas electromagn´eticas en el vac´ıo El problema que se plantea es el siguiente: tenemos una onda electromagn´etica viajando en el vac´ıo que llega a una superficie infinita (la l´amina de grafeno). Queremos averiguar qu´e cantidad de esa onda se refleja y cu´anta se transmite. Para ello, partamos de la base: para describir c´omo se propaga una onda electromagn´etica nos basamos en las ecuaciones de Maxwell (sistema CGS)[11]: ∇·D= 4πρ ∇·B= 0 ∇×E=−1 c ∂B ∂t ∇×H=4π cJ+1 c ∂D ∂t (1) Teniendo en consideraci´on que D=E+ 4πP,H=B−4πMy en el vac´ıo no hay cargas (ρ= 0) ni corrientes (J=0), entonces el sistema se reduce a ∇·E= 0 ∇·H= 0 ∇×E=−1 c ∂H ∂t ∇×H=1 c ∂E ∂t (2) Resolver este sistema es resolver un problema de seis inc´ognitas: las tres componentes del campo el´ectrico Ey las tres del campo magn´etico H. No obstante, se puede comprobar que, en el caso del vac´ıo, el problema se reduce s´olo a dos inc´ognitas: las dos ´ultimas ecuaciones relacionan directamente el rotacional de ambos campos con la derivada temporal del otro, de modo que, si tenemos uno de ellos calculado, el otro puede calcularse directamente a trav´es de esas relaciones. Esto reduce el problema a s´olo tres inc´ognitas (las tres del campo el´ectro o las tres del campo magn´etico). Por otro lado, las dos primeras ecuaciones (divergencia del campo igual a cero) establece una relaci´on directa entre sus tres componentes. Si, por ejemplo, tom´asemos una onda plana1E(r, t) = E0ei(k·r−ωt) tendr´ıamos que2 ∇·E=ˆux ∂ ∂x + ˆuy ∂ ∂y + ˆuz ∂ ∂z E0ei(k·r−ωt)=E0e−iωt(ˆuxikx+ ˆuyiky+ ˆuziky)eik·r= 0, como E0=E0xˆux+E0yˆuy+E0zˆuzy ˆui·ˆuj=δij, entonces se tiene que 1El caso de onda plana es interesante pues, tal como veremos posteriormente, cualquier onda puede expresarse como superposici´on de ondas planas. 2ˆux,ˆuyy ˆuzson los vectores de la base de R3en el espacio eucl´ıdeo. 3 k·E0= 0 ⇒E0z=−kxE0x+kyE0z kz , de modo que, con saber el valor de ky dos componentes del campo el´ectrico, ya conocemos la tercera y, de acuerdo con lo anterior, tambi´en conocemos el campo H. As´ı, nuestro problema consiste en calcular las componentes xeydel campo: tenemos que calcular el vector bidimensional de componentes E0xyE0y. Este vector puede expresarse en una base vectorial de dimensi´on dos; como veremos m´as adelante, la base m´as adecuada es la de modos transversales el´ectricos y magn´eticos, de modo que, teniendo en cuenta que kk=qk2 x+k2 y, definimos los vectores de la base como ap=1 kkkx ky(Modo transversal magn´etico) as=1 kk−ky kx(Modo transversal el´ectrico) (3) Se puede comprobar que, para una onda transversal magn´etica, E0z6= 0 pero que para una transversal el´ectrica E0z= 0. Adem´as, estos vectores componen una base ortonormal3. De esta manera, una onda plana se puede expresar como E=X µ εµaµei(k·r−ωt), donde la suma esta extendida a los dos modos y εµes la componente del campo en cada uno de los vectores de la base. Si ahora aplicamos la expansi´on de ondas planas de Rayleigh [12], podemos expresar cualquier onda como superposici´on de ondas planas. Si cada una de esas ondas tiene una descomposici´on en la base propuesta, una onda cualquiera puede expresarse como E=ZdkX µ εµaµei(k·r−ωt)(4) Una identidad muy ´util que cumplen los modos transversales es la siguiente (probada en el Ap´endice, secci´on 5.1): −ˆuz×Hµ=YµEµ,con Ys=qzeYp=1 qz ,(5) donde hemos definido el vector de ondas normalizado q=k g, con g=qk2 x+k2 y+k2 z=ω c. A las cantidades Yµlas llamamos impedancias. En este caso HµyEµse refieren a los vectores en el plano xy y, adem´as, incluimos la dependencia de onda plana en ellos. Una vez hemos elegido la base en la que trabajaremos y hemos visto sus propiedades, pasamos a estudiar el problema de la transmisi´on y la reflexi´on de ondas a trav´es de una l´amina bidimensional (en nuestro caso, grafeno). 3Es decir, ai·aj=δij . Ver Ap´endice, secci´on 5.1 4 1.2.2. C´alculo de coeficientes de transmisi´on y reflexi´on, confinamiento de modos p La situaci´on que queremos estudiar es c´omo se transmite y refleja una onda que se propaga en el vac´ıo cuando se encuentra con una l´amina bidimensional (el caso del grafeno es este, una capa del grosor de un ´atomo de carbono). Figura 1: Esquema de transmisi´on y reflexi´on de la onda incidente A partir de este punto emplearemos una notaci´on diferente a la habitual, la cual nos facilitar´a notablemente los c´alculos y la lectura. Consideremos lo siguiente: cuando escribimos una onda electromagn´etica de la forma (4) estamos expres´andola en una base, concretamente en la base de ondas planas. As´ı, podemos expresar nuestra onda como la proyecci´on de un estado perteneciente a un espacio de Hilbert Hen la base de ondas planas. Este tratamiento es el mismo que se hace en polarizaci´on a trav´es del c´alculo de Jones [13]. Para cada polarizaci´on tendremos: |k, µi ∈ H 3 hr|k, µi=Eµeik·ryhk, µ|ri=Eµe−ik·r. Adem´as son estados ortonormales: k0, µ0k, µ=k, µ0Zdr|rihr||k, µi=Zdrei(k−k0)Eµ·E0 µ=δµµ0δ(k−k0), donde se ha tomado que I=Zdr|rihr|. De esta manera, la transmisi´on y reflexi´on pueden escribirse en t´erminos de estos elementos del espacio de Hilbert. Puede probarse que cada ket cumple por separado las siguientes relaciones: |E+i=|Eii+|Eri=|k, µ, +i+rµ k|k, µ, −i |E−i=|Eti=tµ k|k, µ, +i (6) de manera que los signos + y - en el ket indican si la onda viaja en sentido positivo o negativo del eje z(es decir, en la dependencia de onda plana tenemos e−ikzzoeikzz). El sistema que referencia que tomamos tiene su origen en la l´amina y toma valores positivos por debajo de ella (por donde viaja la onda transmitida) y valores negativos por encima (desde donde viene la onda incidente y hacia donde va la onda reflejada). Los coeficientes rµ kytµ kson los coeficientes de reflexi´on y transmisi´on respectivamente. Teniendo esto en cuenta y considerando las siguientes ecuaciones de continuidad [11]: 5 ˆuz×(E+−E−) = 0 ˆuz×(H+−H−) = 4π cJ=−4π cσˆuz×(ˆuz×E+), (7) que vienen de la conservaci´on de la componente paralela al plano del campo el´etrico y el salto que tal componente del campo magn´etico sufre debido a la corriente inducida en el plano, puede probarse que, teniendo en cuenta que q=k/g (vector de ondas normalizado en el vac´ıo), α= 2πσ/c (conductividad normalizada) y la relaci´on (5), se llega a las siguientes expresiones para los diferentes modos: rT E,s q=−α α+qz , rT M,p q=−αqz αqz+ 1 tT E,s q=qz α+qz , tT M,p q=1 αqz+ 1 (8) Vamos a probarlo. La primera condici´on no es m´as que la continuidad de la componente paralela a la superficie del campo el´ectrico, lo cual puede expresarse en t´erminos de los coeficientes como 1 + rµ q=tµ q (µ∈ {s, p}). Por otro lado, si partimos de la segunda ecuaci´on: ˆuz×(H+−H−) podemos valernos de la expresi´on −ˆuz×Hµ=YµEµ. No obstante, aqu´ı hay un punto sutil: el valor de la impedancia depende de qu´e signo tenga qz(pues o bien es directamente proporcional a este valor o lo es a su inversa), es decir, de en qu´e sentido viaje la onda en la direcci´on z. El sistema de referencia que nosotros marcamos tiene z= 0 en la placa, por debajo de ella z > 0 y por encima z < 0. De este modo, la ondas ondas incidente y transmitida viajar´an en el sentido positivo del eje z, mientras que la onda reflejada viaja en el sentido negativo (lo cual introduce un signo – en la impedancia). Si tenemos esto en cuenta, poniendo H−=Hi+HryH+=Ht(indicando los ´ındices si es onda incidente i, transmitida to reflejada r): ˆuz×(H+−H−) = ˆuz×(Ht−Hi−Hr) = Yqµ(−Et+Ei−Er). Por otra parte, si usamos la identidad vectorial a×(b×c) = b×(a·c)−c·(a·b), y α=2π cσla conductividad adimensional, se tiene que el otro lado de la ecuaci´on queda como: −4π cσˆuz×(ˆuz×E+) = −2α[ˆuz(ˆuz·E+)−E+] = 2αE+ Donde hemos usado que E+tiene componente znula. Esto se ve en la propia ecuaci´on: si tenemos ese vector igualado a ˆuz×A, con Acualquier vector de R3, el resultado ser´a un vector en el espacio x, y, en R2(pues dar´a un vector mutuamente ortogonal a ˆuzyA). Si juntamos todo lo desarrollado tendremos: Yqµ(−Et+Ei−Er) = 2αE+. Expresemos ahora este resultado en la base de estados |k, µi: Yqµ(−|Eti) + |Eii−|Eri) = 2α|Eti ⇒ Yqµ(−tµ q+ 1 −rµ q)|k, µi= 2αtµ q|k, µi, 6 Donde hemos obviado la parte de + y −del ket, pues se refiere al sentido de propagaci´on de la onda en zy ya lo hemos tenido en cuenta antes. Si proyectamos sobre el bra hk, µ|obtenemos lo siguiente: −tµ q+ 1 −rµ q=2αtµ q Yqµ (9) Como 1 + rµ q=tµ q, si escribimos rµ q=tµ q−1 en (9): −tµ q+ 1 + 1 −tq=2α Yqµtµ q⇒ −tµ q+ 1 = α Yqµtµ q⇒tµ q=Yqµ α+Yqµ Como rµ q=tµ q−1, se tiene que tµ q=Yqµ α+Yqµ rµ q=−α Yqµ+α(10) Si planteamos las impedancias para cada uno de los modos, se obtienen los resultados antes expuestos: rT E,s q=−α α+qz , rT M,p q=−αqz αqz+ 1 tT E,s q=qz α+qz , tT M,p q=1 αqz+ 1 Una de las consecuencias m´as importantes de estos resultados surge al hacerse la siguiente cuesti´on: ¿Es posible tener onda transmitida y reflejada sin que haya onda incidente para alguno de los modos?. De ser as´ı, tales ondas no vendr´ıan de una onda incidente, sino de un plasm´on superficial, una onda que se propaga por la superficie. Si nos planteamos nuevamente las ecuaciones (9) y (7) pero tomando |E+i=rq|k, µ, −i y las relaciones ya calculadas (8) sacamos dos conclusiones: Es imposible que una onda tipo s(transversal el´ectrica) viaje como un plasm´on. La onda tipo p(transversal magn´etica) puede viajar como plasm´on si y s´olo si qz=−1 α. Para sacar la condici´on de plasm´on qz=−1 αbasta con igualar los coeficientes de reflexi´on y transmisi´on (si no hay onda incidente, estos son id´enticos; basta con quitar el 1, que viene de la onda incidente, de la ecuaci´on 1 + rµ q=tµ q). Teniendo todo lo explicado en cuenta, desarrollaremos ahora un modelo unidimensional que trate de explicar c´omo se propaga el plasm´on en la superficie en presencia de defectos variables en el espacio. 7 2. Modelo 2.1. Planteamiento: amplitud como ecuaci´on integral El sistema que se plantea es el siguiente: Figura 2: Esquema del modelo 1D. Planteamos el modelo del siguiente modo: tenemos un campo el´ectrico que viaja a trav´es de la superficie del plasm´on en una direcci´on x. Por simplicidad, plante´emos que se mueve por una red peri´odica unidimensional finita de longitud L, de modo que el m´odulo del campo ser´a de la forma E(x) = eikpx+X G AGei(kp+G)x. Por otro lado, la conductividad (normalizada) ser´a la propia del material (en nuestro caso grafeno) αg4y el aporte que supone el defecto, el cual a˜nade inhomogeneidades en tal conductividad de la forma ∆α(x), siendo este par´ametro la variaci´on relativa de la conductividad. Tales inhomogeneidades pueden ser propias de defectos del material o incluso inducidas a trav´es de un potencial el´ectrico. De este modo, la conductividad del material contando con el defecto ser´a α(x) = αg+ ∆α(x).(11) Si planteamos que J=σE(ley de ohm) y la segunda relaci´on de continuidad de (7), tenemos: −ˆuz×H++ ˆuz×H−= 2α(x)E(x). Considerando que el sistema es sim´etrico respecto al plano que forma la superficie de grafeno5 4La conductividad del grafeno puede obtenerse como σ=σintra +σinter, con σintra =2ie2t ~πΩln 2 cosh 1 2t yσinter = e2 4~1 2+1 πarctan Ω−2 2t−i 2πln (Ω + 2)2 (Ω −2)2+ (2t)2, con Ω = ~ω/µ yt=T/µ, con Ten unidades de energ´ıa. Para m´as referencias consultar [14]. 5Si Exes sim´etrico (que as´ı lo hemos tomado) como la divergencia del campo es nula ∂xEx+∂zEz= 0 ⇒∂zEz=−∂xEx, con lo que Ezes antisim´etrico. El campo, por otro lado, llevar´a el mismo signo que Ez, pues los campos se relacionan con el rotacional. As´ı, el campo Hyser´a positivo por encima de la placa y negativo por debajo. 8 ξ(q) = exp −igq(N−1)γa 2sin gqNγa 2 sin gqγa 2(31) Es importante ver que este factor de estructura ξ(q) surgir´a siempre que tengamos una estructura peri´odica de defectos, de modo que sus propiedades son extrapolables a cualquier tipo de variaci´on en la conductividad siempre que ´esta sea peri´odica. A partir de este resultado podemos estudiar el sistema. El primer paso ser´a considerar la First Order Born Approximation. 3.2. Predicci´on en FOBA Del mismo modo que para una gaussiana, podemos aproximar en FOBA a trav´es de la ecuaci´on (24). De este modo, el coeficiente de reflexi´on en tal aproximaci´on, dado por (25) ser´a R= 2πi α3 gqp B(−qp) 2 ⇒RF OBA =−2πi∆α0(−2qp) α3 gqp ξ(−2qp) 2 . Como |a·b|=|a|·|b|∀a, b ∈C, podemos escribir RF OBA =−2πi∆α0(−2qp) α3 gqp|ξ(−2qp)|2 =−2πi∆α0(−2qp) α3 gqp 2 |ξ(−2qp)|2=RF OBA 0|ξ(−2qp)|2, donde, como hemos expresado en la ecuaci´on, la primera parte ha sido calculada para el caso de una sola gaussiana. Con respecto al t´ermino del factor de estructura: |ξ(−2qp)|2= exp {igqp(N−1)γa}sin (gqpNγa) sin (gqpγa) 2 =|exp {igqp(N−1)γa}| sin (gqpNγa) sin (gqpγa)2 . Como z=|z|eiθ∀z∈C, el m´odulo de un n´umero complejo que consta de una exponencial imaginaria es uno. Por lo tanto, teniendo esto en cuenta, el resutado final ser´a RF OBA =RF OBA 0 sin2(gqpNγa) sin2(gqpγa)(32) 15 A partir de esta primera aproximaci´on podemos comprobar c´omo cambia la reflexi´on en este sistema respecto al de una gaussiana, estudiado en profundidad en ([15]). Si representamos gr´aficamente para N= 1,2,3 y 4 se obtiene el siguiente gr´afico: 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 R a/λp R en FOBA analítica (δ=−0.2,a=20 nm y γ=2) N=1 N=2 N=3 N=4 Figura 5: Coeficiente de reflexi´on en FOBA para N=1,2,3 y 4. En nuestro caso, con N= 1 recuperamos el caso ya conocido. Se aprecia que, adem´as de un pico que crece en magnitud y va haci´endose m´as estrecho conforme Naumenta, tambi´en surgen picos secundarios para longitudes de onda del plasm´on cercanas a la principal. De hecho, a cada lado del pico principal aparecen N−1 m´aximos secundarios. El hecho de que el pico vaya aumentando conforme lo hace Nviene de la indeterminaci´on en el factor de estructura ξ(q) que se da cuando el argumento del seno del denominador se hace nulo. Ve´amoslo: ξ(q) = exp −igq(N−1)γa 2sin gqNγa 2 sin gqγa 2si gqγa 2=mπ, con m∈Z=⇒ξ(q)∝sin(Nmπ) sin(mπ)=0 0. Si hacemos el l´ımite en el cual numerador y denominador tienden a cero: l´ım gqγa−→2mπ sin gqNγa 2 sin gqγa 2≈l´ım gqγa−→2mπ gqNγa 2 gqγa 2 =N, De modo que el factor de reflexi´on Rtiene un pico proporiconal a N2, pues es proporcional al cuadrado del cociente del seno del doble del ´angulo (Que tambi´en converge a Nen la indeterminaci´on). 16 3.3. Resultados de las simulaciones 3.3.1. Discretizaci´on de la ecuaci´on integral Si vamos a la ecuaci´on integral que hemos de resolver para las componentes Fourier B(q) del campo el´ectrico (21), podemos escribirla de este modo: ∆α(q−qp) = −B(q)−Z+∞ −∞ ∆α(q−q0)G(q0)B(q0)dq0,(33) con G(q) = 1 Y(q) + αg (34) la funci´on de green. De este modo, introduciendo la delta de Dirac, que verifica Z+∞ −∞ dxf(x)δ(x− x0) = f(x0), podemos escribir B(q) = Z+∞ −∞ dq0B(q0)δ(q−q0). De este modo, la ecuaci´on (33) queda ∆α(q−qp) = −Z+∞ −∞ dq0B(q0)δ(q−q0)+∆α(q−q0)G(q0)B(q0)=−Z+∞ −∞ dq0B(q0)δ(q−q0)+∆α(q−q0)G(q0). Si ahora discretizamos la ecuaci´on, es decir, Z+∞ −∞ −→ +∞ X q0=−∞ dq0−→ ∆q0 δ(q−q0)−→ δqq0(De delta de Dirac a delta de Kronecker), obtenemos ∆α(q−qp) = − +∞ X q0=−∞ ∆q0B(q0)(δqq0+ ∆α(q−q0)G(q0)).(35) Definimos ahora las siguientes matrices y vectores: Vector F:Fq= ∆(q−qp) Vector B:Bq0=B(q0) Matriz ˜ M:Mqq0= (δqq0+ ∆α(q−q0)G(q0))∆q0. (36) La ecuaci´on (36) se expresa como 17 F=˜ MB(37) De modo que, para cada q(cada iteraci´on), calcularemos la amplitud Ba trav´es de la inversa de ˜ M: B=˜ M−1F.(38) Con esta ecuaci´on ya discretizada podemos programar un c´odigo que la resuelva (ver Ap´endice, apartado 5.8). A la hora de hacerlo, es muy importante determinar c´omo es el integrando. Por ejemplo, la funci´on de Green tiene un polo en qz=p1−q2=−1/αg, de modo que, en tal polo, es necesario disminuir el intervalo de integraci´on, pues la funci´on var´ıa mucho m´as bruscamente. En el caso que nos ocupa (el de Ngaussianas), la ´unica dificultad a˜nadida al integrando es la indeterminaci´on 0 0del factor de estructura para valores de qtales que gqγa = 2mπ. Sabemos que, para esos valores, el factor de estructura vale N(tal y como hemos visto al final de la secci´on 3.2). De ese modo, basta con a˜nadir un condicional en el c´odigo de la siguiente forma a la hora de generar el vector Fy la matriz ˜ M: % F ve ctor and M matrix generator for i =1:2∗N i f rem(gamma∗(real ( q( i ) )−real ( qp ) ) ∗g∗a /2 , pi )==0 F( i ) =1/(4∗sqrt (pi ) ) ∗alphaG∗g∗a∗exp(−(q ( i )−qp ) ˆ2∗( a∗g ) ˆ2/16) ∗delta∗NumGauss ; else F( i ) =1/(4∗sqrt (pi ) ) ∗alphaG∗g∗a∗exp(−(q ( i )−qp ) ˆ2∗( a∗g ) ˆ2/16) ∗delta∗exp(−1 i ∗gamma∗( q (i)−qp) ∗a∗g∗(NumGauss−1)/2) ∗sin (NumGauss∗gamma∗(q( i )−qp ) ∗a∗g /2) / sin (gamma∗(q( i )−qp) ∗a∗g /2) ; end G( i , i )=qz ( i ) /(1+alphaG∗qz ( i ) ) ∗dq( i ) ; G1( i )=qz ( i ) /(1+alphaG∗qz (i)); Q1( i , i ) =1; end % For any reason , M matrix is generated f a s t e r whether i t i s d e fined as the % product of two matrix for i =1:2∗N; for j =1:2∗N; i f rem(gamma∗(real ( q( i ) )−real ( q ( j ) ) ) ∗g∗a /2 , pi )==0 M1( i , j )=−1/(4∗sqrt (pi ) ) ∗alphaG∗a∗g∗exp(−(q ( i )−q( j ) ) ˆ2∗( a∗g ) ˆ2/16) ∗delta∗ NumGauss ; else M1( i , j )=−1/(4∗sqrt (pi ) ) ∗alphaG∗a∗g∗exp(−(q ( i )−q( j ) ) ˆ2∗( a∗g ) ˆ2/16) ∗delta∗exp(−1 i ∗gamma∗(q( i )−q ( j ) ) ∗a∗g∗(NumGauss−1) /2) ∗sin (NumGauss∗gamma∗( q ( i )−q( j ) ) ∗a∗g /2) / sin (gamma∗(q( i )−q ( j ) ) ∗a∗g /2) ; end end end M=M1∗G; } donde rem(a,b) es la funci´on de resto. As´ı, si la cantidad gqγa es un m´ultiplo de π(es decir, rem(gqγa,π)=0), sustituimos el factor de estructura por N(en el caso del c´odigo, NumGauss). 18 3.3.2. Comparaci´on FOBA anal´ıtica y de simulaciones En este apartado comprobaremos que la correcci´on introducida en el c´odigo para el defecto de N gaussianas es correcta: modificamos el c´odigo para que calcule la FOBA, de modo que comentamos todas las operaciones que involucran a la matriz ˜ Me identificamos los vectores F=B, lo cual es equivalente a quitar la integral de la ecuaci´on (notar que ´esta s´olo aparece en la matriz ˜ M). Si comparamos los resultados de las simulaciones con el c´alculo anal´ıtico de la secci´on 3.2, obtenemos lo siguiente: 0 0.01 0.02 0.03 0.04 0.05 0.06 0.07 0.08 0.09 0.1 0 0.1 0.2 0.3 0.4 0.5 0.6 R a/λp R en FOBA analítica (δ=−0.2,a=20 nm y γ=2) comparada con FOBA de simulaciones Simulación N=1 Analítico, N=1 Simulación N=2 Analítico, N=2 (a) N=1,2 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 R a/λp R en FOBA analítica (δ=−0.2,a=20 nm y γ=2) comparada con FOBA de simulaciones Simulación N=3 Analítico, N=3 Simulación N=4 Analítico, N=4 (b) N=3,4 Figura 6: Comparaci´on FOBA anal´ıtica y con simulaci´on Se aprecia una perfecta concordancia de las simulaciones con el resultado anal´ıtico. Con esto comprobamos que la implementaci´on de la modificaci´on para un sistema peri´odico de Ndefectos se ha hecho correctamente en el c´odigo. Pasemos ahora a estudiar casos sin aproximaci´on: inclu´ımos el t´ermino integral en la ecuaci´on. 19 3.3.3. Comparaci´on FOBA y simulaciones A continuaci´on comparamos las predicciones de la FOBA con las simulaciones (incluyendo el t´ermino integral). Los resultados son los siguientes 0 0.005 0.01 0.015 0.02 0.025 0.03 0 0.1 0.2 0.3 0.4 0.5 0.6 R a/λp Comparación con FOBA para a=20 nm y N=1 (δ=−0.2) N=1 simulaciones N=1 FOBA (a) N=1 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.1 0.2 0.3 0.4 0.5 0.6 R a/λp Comparación con FOBA para a=20 nm y N=2 (δ=−0.2) N=2 simulaciones N=2 FOBA (b) N=2 0 0.05 0.1 0.15 0.2 0.25 0 0.1 0.2 0.3 0.4 0.5 0.6 R a/λp Comparación con FOBA para a=20 nm y N=3 (δ=−0.2) N=3 simulaciones N=3 FOBA (c) N=3 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0 0.1 0.2 0.3 0.4 0.5 0.6 R a/λp Comparación con FOBA para a=20 nm y N=4 (δ=−0.2) N=4 simulaciones N=4 FOBA (d) N=4 Figura 7: Comparaci´on de FOBA con los resultados de la simulaci´on Aqu´ı los resultados empeoran: el aspecto de la funci´on es similar, pero se ve un cierto desplazamiento en el eje de abscisas, as´ı como unas formas m´as irregulares en las simulaciones (aunque esto no es debido a la imprecisi´on de la aproximaci´on, sino a la resoluci´on de la simulaci´on). A´un as´ı, se aprecia que la forma funcional es similar (pico resonante de alta magnitud y picos secundarios). Esto nos permite estudiar la aproximaci´on FOBA y traspasar las conclusiones a los casos reales teniendo en cuenta este desplazamiento en las longitudes de onda. 20 3.3.4. Comportamiento con a, anchura del defecto Al igual que en [15], realizamos simulaciones variando la anchura ade los defectos gaussianos. El resultado que obtenemos es el siguiente 0 0.05 0.1 0.15 0.2 0.25 0.3 0.35 0.4 0.45 4 6 8 10 12 14 16 18 20 22 R λ (µm) Reflexión variando la anchura del defecto (N=2;γ=2) a=20 nm; δ= −0.2 a=30 nm; δ= −0.2 a=40 nm; δ= −0.2 a=50 nm; δ= −0.2 a=20 nm; δ= −0.4 a=30 nm; δ= −0.4 a=40 nm; δ= −0.4 a=50 nm; δ= −0.4 (a) Variaci´on con a 0 0.02 0.04 0.06 0.08 0.1 0.12 0 0.1 0.2 0.3 0.4 0.5 0.6 R a/λp R para distintos a (δ=−0.2, N=2, γ=2) a=20 nm a=30 nm a=40 nm a=50 nm (b) Comprobaci´on de ley de escala Figura 8: Variaci´on del coeficiente de reflexi´on Rcon la anchura de las gaussianas a Se aprecia que, conforme aumentamos la anchura a, los picos se desplazan hacia longitudes de onda mayores (resultado que ya se reproduc´ıa en [15] para N=1). Para explicar esto es necesario conocer el origen del pico de reflexi´on: tal pico es debido a la resonancia de la onda que viene del vac´ıo (de longitud de onda λ) con el defecto de anchura a. De este modo, si aumentamos la anchura ala longitud de onda resonante λser´a mayor: las longitudes de onda que resuenan son m´as largas dado que la anchura crece. Por otro lado, si representamos la reflexi´on en funci´on de a/λpse ve que las gr´aficas se superponen entre s´ı. Esto puede apreciarse en la expresi´on de la FOBA: cumple una ley de escala seg´un la magnitud a/λp. Tambi´en se presentan resultados para dos δdiferentes. L´ogicamente, conforme δ(profundidad del defecto) es mayor, la reflexi´on es mayor pues, como vemos en la aproximaci´on FOBA (26), el coeficiente de reflexi´on es proporcional a δ2. Veamos ahora m´as resultados de simulaciones con variaciones en este par´ametro. 21 3.3.5. Comportamiento con δ, altura del defecto Los resultados de los simulaciones para un n´umero de gaussianas N= 2 y 3, anchura del defecto a=20 nm y una distancia relativa entre gaussianas γ= 2 variando δson 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 4 6 8 10 12 14 16 18 20 22 R λ (µm) Reflexión variando δ(a=20 nm; N=2; γ=2) δ=−0.1 δ=−0.2 δ=−0.3 δ=−0.4 δ=−0.7 δ=−0.9 (a) N=2 0 0.2 0.4 0.6 0.8 1 1.2 4 6 8 10 12 14 16 18 20 22 R λ (µm) Reflexión variando δ(a=20 nm; N=3; γ=2) δ=−0.1 δ=−0.2 δ=−0.3 δ=−0.4 δ=−0.7 δ=−0.9 (b) N=3 Figura 9: Variaci´on del coeficiente de reflexi´on Rcon la altura de las gaussianas δ. Vemos que llega un punto de la profundidad (δ=-0,7,-0.9) en el cual la forma de la gr´afica cambia: se va ensanchando el pico central y subiendo su intensidad hasta que rebasa llega justo al 1 (Rno puede ser mayor que 1, pues significa que se refleja m´as energ´ıa de la que entra en el defecto). Un resultado similar se reproduce en [15], donde se usan defectos de anchura mucho mayor (en torno a los micr´ometros). Como veremos en el apartado 4, podemos encontrar aplicaciones en este fen´omeno, las cuales involucran defectos en esa escala. 22 3.3.6. Comportamiento con N, n´umero de gaussianas Las simulaciones para a= 20 nm, δ=-0.2 y γ= 2 para distintos Nson las siguientes: 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 4 5 6 7 8 9 10 R λ (µm) Reflexión variando N(a=20 nm; δ=−0.2;γ=2) N=1 N=5 N=6 N=7 N=8 N=12 Figura 10: Variaci´on del coeficiente de reflexi´on Rcon el n´umero de gaussianas. Los resultados son similares a los que hab´ıamos predicho para FOBA en la parte anal´ıtica: el pico va creciendo conforme aumenta el valor de Ny, adem´as, se va estrechando. No s´olo eso, sino que van apareciendo m´as m´aximos secundarios (N−1 a cada lado), aunque estos var´ıan su intensidad en mucha menor magnitud que el pico central (el correspondiente a la indeterminaci´on 0 0ya comentada). Algo llamativo de estos resultados es que, como ya vimos en 3.2, el m´aximo de la reflexi´on es proporcional aN2, de modo que la reflexi´on puede ser arbitrariamente grande (dado que podemos usar un n´umero N arbitrariamente grande). Esto dar´a lugar a que, para un cierto N, el valor Rrebasar´a la unidad. Podemos estimar a trav´es de la FOBA en qu´e valor de Nocurre esto. Sabemos que, en el m´aximo, el factor de estructura aporta a la reflexi´on un valor N2, de manera que RF OBA m´ax =RF OBA 0,m´ax N2. Si analizamos la funci´on RF OBA 0(ver Ap´endice, secci´on 5.7) podemos extraer que ´esta tiene un m´aximo para a/λp= 1/√2π≈0,22, lo que da lugar a un valor de la reflexi´on RF OBA 0,m´ax ≈0,58δ2. Si elegimos δ=−0,2, tendremos, juntando con el factor de estructura, que RF OBA m´ax (N)=0,0232N2. As´ı, si calculamos este valor para una serie de N: 23 RF OBA m´ax N 0,02 1 0,09 2 0,21 3 0,37 4 0,58 5 0,84 6 1,14 7 1,48 8 1,88 9 2,32 10 Cuadro 1: Valores de RF OBA m´ax em funci´on de N. Como vemos, para el valor de N=7 ya se rebasa la unidad y se pierde el sentido f´ısico del factor de reflexi´on. Sin embargo, nuestras simulaciones llegan hasta N= 12 y no se aprecia esta divergencia. La raz´on es la resoluci´on de la simulaci´on: conforme Naumenta, el pico no s´olo se hace m´as intenso, sino que tambi´en se hace m´as estrecho, de modo que es mucho m´as costoso detectarlo en t´erminos de resoluci´on. De este modo, no lo detectamos porque nuestra discretizaci´on no lo capta, para verlo necesitar´ıamos unas simulaciones con una resoluci´on mucho m´as baja (aqu´ı hemos utilizado δλ ≈0,13nm)8. Por otro lado, podemos dar una visi´on intuitiva de los m´aximos de reflexi´on que tenemos: estos se deben a la resonancia del plasm´on con el defecto, de modo que, si a˜nadimos m´as defectos defectos de la misma anchura (con la cual resuena y provoca el m´aximo), o sea, incrementamos N,la intensidad de esta reflexi´on aumenta, pues resuena con m´as defectos. Por otro lado, la aparici´on de los m´aximos secundarios es una combinaci´on de otros modos de resonancia: el plasm´on puede resonar entre los centros de cada una de las gaussianas, de manera que da lugar a modos mixtos (entre diferentes defectos) de resonancia. 8En realidad la resoluci´on var´ıa seg´un la zona que integremos, tal como coment´abamos al comienzo. De hecho, integramos en la variable adimensional q. 24 con τ=2πi α3 gqpB(qp). Hemos asumido que Imqp>0, de modo que el t´ermino exponencial converger´a. Un tratamiento an´alogo nos llevara a que l´ım x−→−∞ E(x) = eikpx+ρe−ikpx, con ρ=2πi α3 gqp B(−qp). De este modo, los coeficientes de reflexi´on y transmisi´on se definen como R=|ρ|2= 2πi αg B(−qp) 2 T=|1 + τ|2=R=|ρ|2= 1 + 2πi αg B(qp) 2 5.4. C´alculo de integrales con la funci´on de Green Vamos a plantear el c´alculo de una integral del tipo I=Z+∞ −∞ dqG(q)f(q) = Z+∞ −∞ f(q)dq Yq+αg . La funci´on de Green puede expresarse tambi´en como: G(q) = 1 Yq+αg =1 1 qz +αg =qz qz+αg , donde s´olo tenemos en cuenta los modos confinados en la superficie (es decir, los TM). Tenemos que tener en cuenta que qz=p1−q2. Para hacer esta integral nos valdremos del teorema de los residuos de Cauchy, que establece para un recorrido cerrado: If(z)dz= 2πi XRe(f(z0), z0), donde la suma se extiende a los distintos polos situados en z0(distintos puntos) donde la funci´on f(z) tiene polos. Si tenemos en cuenta esto y asumiendo que nuestra f(q) no tendr´a polos, tratemos de calcular Z+∞ −∞ dq 1 + qzαg a trav´es de sus polos. Tal funci´on tiene polos en qzp =−dfrac1αg. A los polos en la variable qlo llamamos qp. Si en la integral hacemos el cambio de variable q=qp−x(As´ı tendremos polos en x= 0: Z+∞ −∞ dq 1 + qzαg =−∞ +∞ dx 1 + p1−(qp−x)2αg, donde hemos empleado que qz=p1−q2. Si asumimos que x << 1, desarrollando la ra´ız: q1−(qp−x)2=q1−q2 p−x2+ 2xqp≈q1−q2 p+ 2xqp. 31 Si ahora sacamos factor com´un (1 −q2 p) = qpz: q1−q2 p+ 2xqp=s(1 −q2 p)(1 −2xqp 1−q2 p ) = qpzs1−2xqp qpz . Considerando que, si x << 1 entocnes √1−x= 1 −x 2, el denominador queda como: 1 + qzp(1 −xqp q2 zp )αg= 1 + qzpαg−xqpαg qzp . Como qzp =−1 αg⇒1 + qzpαg= 0, el numerador queda como −xqpαg qzp y la integral a resolver es Z−∞ +∞ dx −x qzp qpαg =Z+∞ −∞ dx x qzp qpαg . Si le aplicamos el teorema de los residuos en el cero (que es donde est´a el polo): qzp qpαgZ+∞ −∞ dx x=qzp qpαg 2πi l´ım x−→0(x−0)1 x= 2πi qzp qpαg . Si lo juntamos con lo anterior, el resultado es I=Z+∞ −∞ G(q)f(q)dq=2πi qpα3 g f(qp) 5.5. C´alculo de la integral de una gaussiana En este apartado vamos a calcular el valor I=Z+∞ −∞ dxe−x2. Para ello, si consideramos el cuadrado de esta cantidad, tendremos la siguiente integral doble: I2=Z+∞ −∞ dxe−x2Z+∞ −∞ dye−y2. Si pasamos a coordenadas polares, donde r2=x2+y2, dxdy=rdrdθy los extremos de integraci´on pasan a ser rde 0 a +∞yθde 0 a 2π. De este modo: I2=Z+∞ r=0 Z2π θ=0 re−r2drdθ= 2πZ+∞ r=0 re−r2dr= 2π−1 2e−r2+∞ 0 = 2π. O sea que, como I2= 2π, I=Z+∞ −∞ dxe−x2=√2π . 32 5.6. C´alculo de la suma geom´etrica En este apatado vamos a calcular el valor S=PN n=0 rn. Consideremos lo siguiente: S= 1 + r+... +rN rS =r+r2+... +rN+1 Si hacemos ahora rS −S: rS −S=S(r−1) = rN+1 −1⇒S= N X n=0 rn=rN−1−1 r−1 5.7. C´alculo del m´aximo de reflexi´on en FOBA Recordando la forma del coeficiente de reflexi´on en FOBA para una gaussiana: RF OBA O=δ2π 4(kpa)2δ2exp −1 2(kpa)2. Si ahora escribimos kp=2π λy sustituimos luego x=a λp, tenemos: RF OBA 0(λp) = δ2π 42πa λp2 exp (−1 22πa λp2)⇒RF OBA 0(x) = δ2π3xexp {−2πx}. Si ahora derivamos respecto de xe igualamos a cero: ∂RF OBA 0(x) ∂x = 0 ⇒π3δ3(e−2πx −2πxe−2πx) = 0 ⇒1−2πx = 0 ⇒x=1 2π. Si deshacemos el cambio de variable: x=1 2π⇒a λp =1 √2π. Aqu´ı tenemos un extremo. Para saber si es un m´aximo o un m´ınimo deber´ıamos hacer la segunda derivada. No obstante, si sustituimos en la expresi´on y ´esta no se anula (pues el valor m´ınimo que alcanza el factor de reflexi´on es cero), ser´a un m´aximo: RF OBA 0a λp =1 √2π=δ2π2 2 1 e≈0,58δ2. 33 5.8. C´odigo MatLab Aqu´ı se presenta el c´odigo MatLab empleado para las simulaciones, concretamente para el caso N= 2, a = 20 nm, δ =−0,2 y γ= 2. Este c´odigo ha sido desarrollado (para el caso de un defecto gaussiano) por Pablo Pons9, estudiante de doctorado en el departamento de F´ısica de la Materia Condensada en la Universidad de Zaragoza. Las modificaciones realizadas para el caso de las Ngaussianas se han comentado en la secci´on (3.3.1).Las simulaciones se han realizado en el cluster del Instituto de investigaci´on de biocomputaci´on y F´ısica de Sistemas Complejos (BiFi) con distintos c´odigos como ´este, variando los distintos par´ametros que se han estudiado en la secci´on 3.3. function main %= echo ( del t a , a ) % Phy s ical consta nts c =2.99792458 e8 ; %[m/s ] speed of l i g h t h=4.135667516 e −15; %[ eV ] planck ’ s co ns tan t eps0 =8.8541878176 e −12; %[F/m] vacuum p e r m i t i v i t t y e=−1.602176565e−19; %[C] e l e c t r o n charge % %FEM Parameters N1=1000; % number o f v al ue s c a l c u l a t e d between qmax and qp+e1 N2=600; % . . . qp+e2 and qp−e2 ; should be an even number to s p l i t p oint s in two regions N22=600; % . . . N3=300; % . . . qp−e2 and 1+e1 N4=300; % . . . 1+e1 and 1−e1 ; should be an even number to avoid 1 zero N5=300; % . . . 1−e1 and 0 N=N1+N2+N22+N3+N4+N5 ; % e1 =0.33; % h a l f range f or zone 4 e2 =4; % h a l f range f or zone 2 e3=150; % times the r ea l part of the pole f or h a l f range f or zone 2.2 f3=1000000; % qmf=1; % % % Graphene parameters mu=0.2; %[ eV] chemical p o t e n t i a l del t a =−0.2; %d e l t a=str2num ( d e l t a ) ; % d e f e c t Depth [ −1:0] a=20e −9; %a=str2num ( a ) ; %[m] d e f e c t Width % %Many gaussian parameters NumGauss=2; %Number of gaussian d i s t r i b u t i o n s we are using . gamma=2; %The d is ta nc e between two a djacent gaussia ns i s d=gamma∗a , thus gamma measure how much are each gaussian sep arated one from another . % Spectra range Nk=130; % number of p oint s ev a l u a ted ( at l e a s t 2) lmin=5e −6; % minimum vacuum wavelength lmax=22e −6; % maximum vacuum wavelength % % Relexion and vacuum wavelength vec t o r d e f i n i t i o n R=zeros(Nk : 1 ) ; % R e f l e xi on ( should be [ 0 : 1 ] ) 9Contacto: p[email protected] 34 l f=zeros(Nk: 1 ) ; % Vacuum wavelength % sout=sprintf( ’C:/TFG/d0. %da %dN%dg %d R . txt ’ ,−delta∗10 , a ∗10ˆ9 ,NumGauss ,gamma) ; sbq=sprintf( ’C:/TFG/d0. %da %dN%dg %d Bq . txt ’ ,−delta∗10 , a ∗10ˆ9 ,NumGauss ,gamma) ; % Opening f i l e s output = fopen( sout , ’w ’ ) ; Bq =fopen( sbq , ’w ’ ) ; % for k=1:Nk % Operation Point lambda=lmin+(lmax−lmin ) /(Nk−1) ∗(k−1) ; % vacuum wavelength omega=c/lambda∗2∗pi ;% angular frequency g=omega/c ; % vacuum wavevector sigma=pi ∗e ˆ2/(2∗h∗abs ( e ) ) ∗sqrt (−1) ∗(8∗mu/(h∗omega) −1/(2∗pi )∗log ((2∗mu+h∗omega /(2∗pi ) ) ˆ2/(2∗mu−h∗omega /(2∗pi ) ) ˆ2) ) ; %c o n d u c t i v i t y alphaG=sigma /(2∗eps0 ∗c ) ; % normalized c o n d u c t i v i t y alphaG=1 i ∗imag( alphaG )+imag( alphaG )/ f3 ; % qp=sqrt(1−1/alphaG ˆ2) ; % graphene normalized wavevector ( withou t d e f e c t s ) qpz=imag(sqrt(1−qp ˆ2) ) /abs (imag(sqrt(1−qpˆ2) ) ) ∗sqrt(1−qpˆ2) ; % qpmax=qmf∗20∗sqrt (2 ) /( a∗g )+2∗real ( qp ) ; % maximum wavevector % % Matrix/ ve c tor d e f i n i t i o n M1=zeros(2∗N, 2∗N) ; % M=zeros(2∗N,2∗N) ; %M matrix (M1∗G) F=zeros(2∗N, 1 ) ; % Independent terms ve c t o r (FOBA) B=zeros(2∗N, 1 ) ; % Fi e l d v e c t or G=zeros(2∗N,2∗N) ; % G1=zeros(2∗N, 1 ) ; % Q1=zeros(2∗N,2∗N) ; % q=zeros(2∗N, 1 ) ; % normalized wavevector qz=zeros(2∗N, 1 ) ; % dq=zeros(2∗N, 1 ) ; % % % q and dq v ec to r s generator for i =1:N1 % zone 1 q( i )=−qpmax+(qpmax−real (qp)−e2 ) /N1∗( i −0.5) ; q( i+N+N5+N4+N3+N2+N22)=real (qp )+e2+(qpmax−real ( qp )−e2 ) /N1∗( i −0.5) ; dq( i )=(qpmax−real (qp )−e2 ) /N1 ; dq( i+N+N5+N4+N3+N2+N22)=(qpmax−real (qp)−e2 ) /N1 ; end for i =1:N2/2 % zone 2 q( i+N1)=−real ( qp )−e2+(e2−abs (imag( qp ) ) ∗e3 ) /N2∗2∗( i −0.5) ; q( i+N1+N2/2+N22)=−real ( qp)+abs (imag( qp ) ) ∗e3+(e2−abs (imag( qp ) ) ∗e3 ) /N2∗2∗( i −0.5) ; q( i+N+N5+N4+N3)=real ( qp )−e2+(e2−abs (imag( qp ) ) ∗e3 ) /N2∗2∗( i −0.5) ; q( i+N+N5+N4+N3+N2/2+N22)=real (qp )+abs (imag( qp ) ) ∗e3+(e2−abs (imag( qp ) ) ∗e3 ) /N2∗2∗( i −0.5) ; dq( i+N1)=(e2−abs (imag( qp ) ) ∗e3 ) /N2∗2; dq( i+N1+N2/2+N22)=(e2−abs (imag( qp ) ) ∗e3 ) /N2∗2; dq( i+N+N5+N4+N3)=(e2−abs (imag( qp ) ) ∗e3 ) /N2∗2; dq( i+N+N5+N4+N3+N2/2+N22)=(e2−abs (imag( qp ) ) ∗e3 ) /N2∗2; end for i =1:N22 % zone 22 q( i+N1+N2/2)=−real ( qp)−abs (imag( qp) ) ∗e3+2∗abs (imag(qp) ) ∗e3 /N22∗( i −0.5) ; q( i+N+N5+N4+N3+N2/2)=real ( qp )−abs (imag( qp ) ) ∗e3+2∗abs (imag( qp ) ) ∗e3/N22∗( i −0.5) ; dq( i+N1+N2/2)=2∗abs (imag( qp ) ) ∗e3 /N22 ; dq( i+N+N5+N4+N3+N2/2)=2∗abs (imag( qp ) ) ∗e3 /N22 ; end for i =1:N3 % zone 3 35 q( i+N1+N2+N22)=−real ( qp )+e2+(real ( qp )−e2−1−e1 ) /N3∗( i −0.5) ; q( i+N+N5+N4)=1+e1+(real ( qp )−e2−1−e1 ) /N3∗( i −0.5) ; dq( i+N1+N2+N22)=(real ( qp )−e2−1−e1 ) /N3 ; dq( i+N+N5+N4)=(real ( qp)−e2−1−e1 ) /N3; end for i =1:N4 % zone 4 q( i+N1+N2+N22+N3)=−1−e1+2∗e1 /N4∗( i −0.5) ; q( i+N+N5)=1−e1+2∗e1 /N4∗( i −0.5) ; dq( i+N1+N2+N22+N3)=2∗e1 /N4; dq( i+N+N5)=2∗e1/N4 ; end for i =1:N5 % zone 5 q( i+N1+N2+N22+N3+N4)=−1+e1+(1−e1 ) /N5∗( i −0.5) ; q( i+N)=(1−e1 ) /N5∗( i −0.5) ; dq( i+N1+N2+N22+N3+N4)=(1−e1 ) /N5 ; dq( i+N)=(1−e1 ) /N5 ; end % %qz ve c tor generator for i =1:2∗N qz ( i )=sqrt(1−q ( i ) ˆ2) ; end % %Theory % % Integral % B( q )=−\Delta\alpha ( q−q p )−\ int {−\ infty}ˆ{+\infty}\ Delta\alpha (q−q ’ )G( q ’ )B( q ’ ) dq ’ % %FEM %\Delta\alpha (q−q p )=\sum {q’=−q{max}}ˆ{q{max}}[−\ delta {q , q ’}−\ Delta\alpha ( q−q ’ )G( q ’ ) \Delta q ’ ] B(q ’ ) % F=[M]∗B % % F ve ctor and M matrix generator for i =1:2∗N i f rem(gamma∗(real ( q( i ) )−real ( qp ) ) ∗g∗a /2 , pi )==0 F( i ) =1/(4∗sqrt (pi ) ) ∗alphaG∗g∗a∗exp(−(q ( i )−qp ) ˆ2∗( a∗g ) ˆ2/16) ∗delta∗NumGauss ; else F( i ) =1/(4∗sqrt (pi ) ) ∗alphaG∗g∗a∗exp(−(q ( i )−qp ) ˆ2∗( a∗g ) ˆ2/16) ∗delta∗exp(−1 i ∗gamma∗( q (i)−qp) ∗a∗g∗(NumGauss−1)/2) ∗sin (NumGauss∗gamma∗(q( i )−qp ) ∗a∗g /2) / sin (gamma∗(q( i )−qp) ∗a∗g /2) ; end G( i , i )=qz ( i ) /(1+alphaG∗qz ( i ) ) ∗dq( i ) ; G1( i )=qz ( i ) /(1+alphaG∗qz (i)); Q1( i , i ) =1; end % For any reason , M matrix i s generated f a s t e r whether i t i s d efin e d as the % product of two matrix for i =1:2∗N; for j =1:2∗N; i f rem(gamma∗(real ( q( i ) )−real ( q ( j ) ) ) ∗g∗a /2 , pi )==0 M1( i , j )=−1/(4∗sqrt (pi ) ) ∗alphaG∗a∗g∗exp(−(q ( i )−q( j ) ) ˆ2∗( a∗g ) ˆ2/16) ∗delta∗ NumGauss ; else M1( i , j )=−1/(4∗sqrt (pi ) ) ∗alphaG∗a∗g∗exp(−(q ( i )−q( j ) ) ˆ2∗( a∗g ) ˆ2/16) ∗delta∗exp(−1 i 36 ∗gamma∗(q( i )−q ( j ) ) ∗a∗g∗(NumGauss−1) /2) ∗sin (NumGauss∗gamma∗( q ( i )−q( j ) ) ∗a∗g /2) / sin (gamma∗(q( i )−q ( j ) ) ∗a∗g /2) ; end end end M=M1∗G; for i =1:2∗N M( i , i )=M( i , i ) −1; end % %Solve f i e l d B % B=M\F; % % S1=0; for i =1:N4+N5 S1=S1+dq (N−(N4+N5)/2+ i ) /qz (N−(N4+N5)/2+ i ) ∗abs (G1(N−(N4+N5)/2+ i ) ∗B(N−(N4+N5)/2+ i ) ) ˆ2; end S1=S1∗4∗pi/real (qp) /abs ( alphaG ) ˆ3 % % l f (k )=lambda ; B(N1+N2/2+N22/2) B(2∗N−N1−N2/2+N22/2+1) R1=abs(−2∗pi ∗(B(N1+N2/2+N22/2)+B(N1+N2/2+N22/2+1) ) /2∗sqrt (−1) /( alphaG ˆ3∗qp) ) ˆ2 T1=abs(1−2∗pi ∗(B(2∗N−N1−N2/2−N22/2)+B(2∗N−N1−N2/2−N22/2−1) ) ∗sqrt (−1) /2/( alphaGˆ3∗qp ) ) ˆ2 R(k )=R1 ; plot (real (q ) , real (B) , real ( q) ,imag(B) , real ( q ) ,imag(F) ) % % Writing f i l e s for i =1:2∗N fprintf(Bq, ’ %e\t %e\t %e\t %e\t %e\n ’ , l f ( k) , q ( i ) , real (B( i ) ) , imag(B( i ) ) , imag(F( i ) ) ) ; end fprintf(Bq , ’ \n ’ ) ; fprintf( output , ’ %e\t %e\t %e\t %e\t %e\n ’ , l f (k ) , real ( qp ) /g , R( k ) , T1 , S1 ) ; % end plot ( l f ,R) disp ( ’ Simulation i s over ! ’ ) ; % Closing f i l e s fclose( output ) ; fclose(Bq) ; % end 37 Referencias [1] Debije, P; Scherrer, P. Interferenz an regellos orientierten Teilchen im R¨ontgenlicht I. Physikalische Zeitschrift 1916, 17: 277. [2] Wallace, P. R.The Band Structure of Graphite. Physical Review 1947, 71: 622?634 [3] Mouras, S.; et al. Synthesis of first stage graphite intercalation compounds with fluorides.Revue de Chimie Minerale 1987, 24: 572. [4] Ruess, G.; Vogt, F. H¨ochstlamellarer Kohlenstoff aus Graphitoxyhydroxyd. Monatshefte f¨ur Chemie 1948, 78 (3 4): 222. [5] Boehm, H. P.; Clauss, A.; Fischer, G.; Hofmann, U. Proceedings of the Fifth Conference on Carbon. Pergamon Press. 1962 [6] Oshima, C.; Nagashima, A. Ultra-thin epitaxial films of graphite and hexagonal boron nitride on solid surfaces. J. Phys.: Condens. Matter 1997. [7] Novoselov, K. S.; Geim, A. K.; Morozov, S. V.; Jiang, D.; Zhang, Y.; Dubonos, S. V.; Grigorieva, I. V.; Firsov, A. A. Electric Field Effect in Atomically Thin Carbon Films Science 2004, 306 (5696): 666?669. arXiv:cond-mat/0410550. Bibcode:2004Sci...306..666N. doi:10.1126/science.1102896. PMID 15499015. [8] Tsoukleri G.; Parthenios, J., Papagelis,K.; Jalil, R., Ferrari, A.C.;Geim,A.K; Novoselov K. S. and Galiotis, C. Subjecting a graphene monolayer to tension and compression. [9] Castro Neto et al.The electronic properties of graphene [10] Pop et al.,Thermal properties of graphene: Fundamentals and applications. [11] Jackson J.D. Classical electrodynamics. 3aedici´on. EEUU: John Wiley & Sons, 1962. [12] S´anchez-Gil, J.A.; Maradudin, A.A. Near-Field and Far-Field Scattering of Surface Plasmon Polaritons by One-Dimensional Surface Defects. Phys. Rev. B 1999, 60, 8359 8367. [13] E. Collett, Field Guide to Polarization, SPIE Field Guides vol. FG05, SPIE (2005). ISBN 0-81945868-6. [14] Nikitin, A. Yu.;Guinea, F.; Garc´ıa-Vidal, F.J.; Mart´ın-Moreno, L., Fields radiated by a nanoemitter in a graphene sheet. Physical Review B 84, 195446, 2011, App A. [15] Garc´ıa-Pomar, J.L; Yu. Nikitin A y Martin-Moreno L. Scattering of Graphene Plasmons by Defects in the graphene Sheet. 2013. [16] Nagase, M; Hibino, H; Kageshima, H.; Yamaguchi, H. Local conductance measurements of DoubleLayer Graphene on SiC substrate. Nanotechnology 2009, 20, 445704 [17] Gorbachev, R.V.; Mayorov, A. S.; Savchenko, A.K.; Horsell, D.W.; Guinea, F. Conductance of p-n-p Graphene Structures with Air-Bridge Top Gates. Nano Lett. 2008, 8, 1995 1999. [18] Liu, G; Jairo Velasco, J.; Bao, W.; Lau, C.N. Fabrication of Graphene p-n-p Junctions with Contactless Top Gates Appl. Phys. Lett. 2008,92,203103. 38