Repositorio Institucional de Documentos

Un problema de bosones que interaccionan. C´alculos anal´ıticos y num´ericos. Jorge Calvo Ibar Dirigido por: David Zueco y Uta Naether 24 de junio de 2014 ii ´ Indice general 1. Introducci´on 1 1.1. Objetivos .................................................... 2 2. Oscilador Arm´onico 3 2.1. OsciladorArm´onicoCl´asico.......................................... 3 2.2. OsciladorArm´onicoCu´antico......................................... 4 3. Modelos Tight Binding Bos´onicos 7 3.1. El m´etodo de c´alculo de bandas Tight Binding - Funciones de Wannier. . . . . . . . . . . . . . . . . . 7 3.2. Bandas de energ´ıa con el modelo Tight Binding para el caso de cadena infinita de osciladores acoplados (1D)........................................................ 8 4. Modelos Tight Binding en una topolog´ıa Diente de Sierra 11 4.1. Diente de sierra, caso de frecuencias iguales. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.2. Diente de sierra, caso de frecuencias distintas. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13 4.3. C´alculo de las bandas planas para dos frecuencias distintas . . . . . . . . . . . . . . . . . . . . . . . . 15 4.4. C´alculo de las bandas planas para casos m´as particulares . . . . . . . . . . . . . . . . . . . . . . . . . 20 5. Diente de sierra con interacci´on no lineal 23 5.1. SingleExcitation................................................ 23 5.2. M´etodonum´erico1............................................... 27 5.3. M´etodonum´erico2............................................... 28 5.3.1. Construyendolabase: ......................................... 28 5.3.2. Construyendolamatriz:........................................ 29 5.4. Comprobaci´on de los m´etodos num´ericos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 5.5. Resultados.................................................... 34 5.5.1. Casolineal ............................................... 34 5.5.2. Casonolineal.............................................. 37 6. Conclusiones 41 A. C´alculo de bandas de cadena infinita de osciladores acoplados 43 B. C´alculo de las bandas en topolog´ıa de tipo dientes de sierra para frecuencias iguales 45 C. C´alculo de las bandas en topolog´ıa de tipo dientes de sierra para frecuencias distintas 49 D. M´etodo Alternativo de c´alculo de la condici´on de Banda Plana 53 E. M´etodo num´erico 1 55 E.1.Explicaci´ondelprograma ........................................... 55 E.2.Resumendelc´odigo .............................................. 58 F. M´etodo num´erico 2 59 F.1.Explicaci´ondelprograma ........................................... 59 F.2.Resumendelc´odigo .............................................. 63 G. C´odigo M´etodo 2 - Modificaci´on para comprobar los auto estados localizados 65 iii H. C´odigo M´etodo 2 - Modificaci´on para comprobar los auto estados localizados, frecuencias distintas 67 I. C´odigo M´etodo 2 - Modificaci´on para el caso no lineal 69 iv ´ Indice de figuras 3.1. Potencial peri´odico al que se someten las part´ıculas y sus funciones de onda de part´ıculas localizadas. (Imagenmodificadade[4]) .......................................... 7 3.2. Cadenadeosciladores.............................................. 8 3.3. Bandas de la cadena de osciladores acoplados. (Notar que k∈[−π, π], la primera zona de Brillouin.) . 9 4.1. Disposici´on en diente de sierra de los sitios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11 4.2. Bandas de energ´ıa Ω = 2, J0= 1, J= 2 y ωa= 3 con +→azul y −→rojo ............. 14 4.3. Bandas de energ´ıa J= 5,∆Ω = −2y ωa= 6 con +→azul y −→rojo ............... 17 4.4. Bandas de energ´ıa J=−5,∆Ω = −2y ωa= 6 con +→azul y −→rojo .............. 18 4.5. Bandas de energ´ıa J∈[−10,10],∆Ω = −2y ωa= 6 con +→rojo y −→azul........... 19 4.6. Bandas de energ´ıa J∈[−10,10], k ∈[−π, π], ω = 1 y a = 0,1 con +→rojo y −→azul . . . . . 20 4.7. Bandas de energ´ıa J∈[−10,10], k ∈[−π, π], ω = 1 y a = 2 con +→rojo y −→azul . . . . . . 20 5.1. Disposici´on en diente de sierra de los sitios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23 5.2. Bandas de energ´ıa J= 5,∆Ω = −2y ωa= 6 con +→azul y −→rojo ............... 25 5.3. Bandas de energ´ıa obtenidas num´ericamente J= 5,∆Ω = −2y ωa= 6 con L´ınea Verde −→ Nivel de energ´ıa m´aximo de la banda de energ´ıa superior, L´ınea Rosa −→ Nivel de energ´ıa m´ınimo de la banda de energ´ıa superior, L´ıneas Azul y Rojo −→ niveles de energ´ıa intermedios que aparecen debido a condiciones de frontera libres y L´ınea Negra −→ Nivel de energ´ıa de la banda plana. . . . . . . . . . 26 5.4. Gap de energ´ıa J= 5,∆Ω = −2y ωa=6.................................. 26 5.5. Unautoestadolocalizado. ........................................... 34 5.6. Autovalores para el caso de 7 sitios y 2 excitaciones calculados por el M´etodo 2 en funci´on del par´ametro no lineal “U”. Para valores de ωa= 6, J=−5 y ∆Ω = −2. Aplicando las condiciones para J0para que existaunabandaplana. ............................................ 37 5.7. Autovalores para el caso de 7 sitios y 2 excitaciones calculados por el M´etodo 2 en funci´on del par´ametro no lineal “U”. Para valores de ωa= 6, J=−5 y ∆Ω = −2. Aplicando las condiciones para J0para que existaunabandaplana. ............................................ 38 5.8. Autovalores para el caso de 7 sitios y 3 excitaciones calculados por el M´etodo 2 en funci´on del par´ametro no lineal “U”. Para valores de ωa= 6, J=−5 y ∆Ω = −2. Aplicando las condiciones para J0para que existaunabandaplana. ............................................ 38 5.9. Autovalores para el caso de 9 sitios y 2 excitaciones calculados por el M´etodo 2 en funci´on del par´ametro no lineal “U”. Para valores de ωa= 6, J=−5 y ∆Ω = −2. Aplicando las condiciones para J0para que existaunabandaplana. ............................................ 39 5.10. Autovalores para el caso de 9 sitios y 3 excitaciones calculados por el M´etodo 2 en funci´on del par´ametro no lineal “U”. Para valores de ωa= 6, J=−5 y ∆Ω = −2. Aplicando las condiciones para J0para que existaunabandaplana. ............................................ 39 A.1.Cadenadeosciladores.............................................. 43 B.1. Disposici´on en diente de sierra de los sitios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45 C.1. Disposici´on en diente de sierra de los sitios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 49 E.1. Disposici´on en diente de sierra de los sitios. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 55 v vi Cap´ıtulo 1 Introducci´on En f´ısica aprendemos a modelizar de la forma m´as simple posible los fen´omenos naturales. Intentamos capturar lo esencial y escribirlo con letras griegas en un Hamiltoniano. Sabemos que nuestro modelo es aproximado y efectivo, pero nos damos por satisfechos si explica el experimento. Pero, por muy simples que sean, los modelos relevantes rara vez son f´aciles de resolver. M´as all´a del oscilador arm´onico o un sistema de dos estados pocos problemas admiten una soluci´on exacta. Es por ello que la f´ısica, adem´as de esa b´usqueda de la simplicidad, sea tambi´en una b´usqueda de aproximaciones. Teor´ıa de perturbaciones, aproximaciones arm´onicas, l´ıquidos de Fermi, etc ... conforman una serie de herramientas que hacen resolubles esos modelos simples, pero dif´ıciles a la vez. Una alternativa a las aproximaciones es la resoluci´on num´erica. Simulaciones por ordenador permiten explorar modelos m´as all´a de reg´ımenes perturbativos as´ı como comprobar las aproximaciones. La combinaci´on de simulaci´on y aproximaciones forman la vida de la mayor´ıa de los f´ısicos te´oricos. Como la vida no puede ser tan sencilla, las simulaciones no son la soluci´on final. Problemas muy importantes son demasiado “grandes” para solucionarlos con ordenador. Cuando decimos “grandes” queremos decir que el n´umero de grados de libertad necesarios son tantos que los ordenadores no pueden con ellos. Un ejemplo t´ıpico es el plegamiento de prote´ınas. En mec´anica cu´antica este problema es todav´ıa m´as grave. En materia condensada, por ejemplo, estamos acostumbrados a trabajar con objetos con muchos grados de libertad (ej. un s´olido). El ´exito de esta f´ısica, inherentemente many body, ha sido desarrollar teor´ıas donde los electrones son esencialmente libres. Si queremos atacar problemas donde las interacciones juegan un papel estamos perdidos. El problema crece exponencialmente con el n´umero de grados de libertad, y ¡tenemos muchos! Para explicar esto un poco mejor, debemos pensar en que cada, por ejemplo, electr´on “vive” en un espacio de Hilbert de dimensi´on D. Para tratar Nelectrones usamos el producto tensorial que sabemos crece de manera exponencial como DN. Es por ello, que dichos problemas no se puedan simular con ordenador [1]. La visi´on anterior es muy negativa. En realidad la din´amica de esos electrones interactuantes no va a visitar ese enorme espacio de Hilbert. T´ıpicamente solo visitar´a una parte exponencialemente peque˜na del mismo. Es por eso, que no debe obsesionarnos el tama˜no sino identificar ese peque˜no subespacio relevante. Debemos buscar simetr´ıas, cantidades conservadas, f´ısica esperable, etc ... Si identificamos la parte relevante somos capaces de resolver el problema [2]. En este Trabajo Fin de Grado presentamos un acercamiento num´erico a un problema dif´ıcil. Se trata de bosones interactuantes que viven en una topolog´ıa entre una y dos dimensiones. El problema tiene su inter´es en problemas de frustraci´on. Explicamos las cantidades conservadas del problema y generamos varios c´odigos que permiten resolver el modelo exactamente en redes no tan peque˜nas. Discutimos la teor´ıa de bandas del modelo en el r´egimen lineal, donde es resoluble anal´ıticamente. Encontramos los par´ametros del modelo donde existen bandas planas. Estudiamos el papel de la interacci´on en estas bandas planas. 1 1.1. Objetivos Los principales objetivos que abordaremos a lo largo de este trabajo son los siguientes: 1. Estudiaremos un sistema de bosones interactuantes tanto el caso lineal como el no lineal. 2. En el caso lineal obtendremos las bandas de energ´ıa y las condiciones que han de cumplirse para su existencia. 3. Estudiaremos el sistema para diferente n´umero de tama˜nos del mismo as´ı como de excitaciones. 4. La resoluci´on del caso no lineal nos obligar´a a obtener un m´etodo num´erico eficiente que permita obtener la dependencia del sistema con la no linealidad. 5. Todo esto permitir´a obtener propiedades del sistema como la degeneraci´on de los autoestados localizados de las bandas planas, la dependencia con el tama˜no del sistema etc. 2 Cap´ıtulo 2 Oscilador Arm´onico [3] En este trabajo centraremos principalmente nuestra atenci´on en un modelo de bosones acoplados y sobre todo la resoluci´on num´erica del mismo. Por completitud estudiaremos a modo introductorio, el caso de un solo bos´on, para as´ı tener una visi´on general del problema el sistema m´as conocido en f´ısica, el oscilador arm´onico, empezaremos brevemente recordando el caso cl´asico para a continuaci´on meternos directamente el caso cu´antico. 2.1. Oscilador Arm´onico Cl´asico Un oscilador arm´onico cl´asico es aquel cuerpo cuyo movimiento viene determinado por el ya conocido potencial cuadr´atico V(x), funci´on del la posici´on xdel mismo1: V(x) = 1 2kx2(2.1) Por lo tanto usando la relaci´on entre fuerza y potencial y la segunda ley de Newton podemos obtener la ecuaci´on diferencial que determinar´a el movimiento: md2x dt2=−dV dx =−kx (2.2) Resolviendo la ecuaci´on diferencial de coeficientes constantes se obtiene la ecuaci´on del movimiento: x=xMcos(ωt −φ) con ω=rk m(2.3) En donde xMes la amplitud del movimiento y ωsu frecuencia. Vamos a obtener su energ´ıa, sumando energ´ıa cin´etica y potencial y observaremos que no depende del tiempo puesto que nos encontramos en un sistema conservativo: E=T+V=1 2mdx dt 2 +1 2mω2x2=··· =1 2mω2x2 M(2.4) Que como se puede ver al sustituir la soluci´on de la ecuaci´on diferencial se obtiene una energ´ıa que se mantiene constante en el tiempo tal y como es propia de los sistemas conservativos. Si ahora fijamos un valor de la energ´ıa total Ey expresamos la energ´ıa total en funci´on del momento py la posici´on xpodemos escribir: E=p2 2m+1 2mω2x2=p2 a2+x2 b2con: a=√2m, b =r2 mω2(2.5) Es decir en el espacio de pyxes una elipse de semiejes aybde manera que cuando la posici´on es m´axima (xM) el momento es m´ınimo (p= 0), o lo que es lo mismo cuando la energ´ıa potencial es m´axima la cin´etica es m´ınima e igualmente a la inversa. 1Consideraremos el caso del oscilador unidimensional para hacer m´as simple la discusi´on. 3 10 Cap´ıtulo 4 Modelos Tight Binding en una topolog´ıa Diente de Sierra En este cap´ıtulo nos centraremos en la obtenci´on te´orica de las bandas de energ´ıa en un modelo tipo Tight Binding con una distribuci´on de osciladores en forma de diente de sierra. Obtendremos las bandas de manera te´orica de un modelo de osciladores acoplados linealmente, para luego desarrollar un m´etodo num´erico que permita tambi´en obtenerlo. Este desarrollo nos permitir´a luego poder incluir interacciones no lineales que no son accesibles para una resoluci´on anal´ıtica. Los c´alculos te´oricos del modelo lineal nos permitir´an compararlo con el num´erico y comprobar el funcionamiento correcto del m´etodo num´erico planteado. Hemos escogido una disposici´on en diente de sierra porque es un arreglo a caballo entre una disposici´on en una dimensi´on, como hemos visto en la secci´on anterior, y una red geom´etrica en dos dimensiones; adem´as de contener propiedades interesantes como bandas planas. En la literatura tambi´en es conocido por mostrar el fen´omeno de la frustraci´on del esp´ın. 4.1. Diente de sierra, caso de frecuencias iguales. Consideremos tal y como hemos dicho que tenemos nuestros osciladores en una disposici´on de diente de sierra. En este caso llamaremos a los operadores creaci´on/destrucci´on de la fila superior como a† jyajy a los de la fila de abajo como b† jybj. Consideraremos adem´as que la frecuencia de oscilaci´on de los osciladores de la fila superior y la inferior es la misma, ωya trataremos el caso de frecuencias distintas en la siguiente secci´on. En cuanto al acoplamiento, ser´a de car´acter lineal igual que el usado en la secci´on anterior con dos constantes de acoplamiento: J0que regula el acoplamiento entre los osciladores de la fila superior con la inferior y Jque regula el acoplamiento entre los osciladores de la fila inferior, como se muestra en la siguiente figura: Figura 4.1: Disposici´on en diente de sierra de los sitios. Con todo lo dicho hasta ahora el Hamiltoniano a trata, del cual obtendremos las bandas de energ´ıa: H= N−1 X j=0 ω(a† jaj+b† jbj)+J N−1 X j=0 b† jbj+1 +b† j+1bj+J0 N−1 X j=0 a† jbj+a† jbj+1 +b† jaj+b† j+1aj(4.1) Notar que el car´acter lineal de los acoplamientos permite que se siga conservando que el operador Hamiltoniano conmute con el operador n´umero Nj=a† jaj(igual para las b0s) cumpli´endose las relaciones de conmutaci´on (3.5) as´ı como las combinaciones de operadores a0syb0s([aj, bj0] = [a† j, bj0]=[aj, b† j0] = 0, as´ı como en el espacio Fourier, con demostraciones an´alogas a la secci´on anterior). 11 As´ı al tratarse de un modelos con acoplamientos a primeros vecinos tipo Tight Binding o de enlace fuerte podemos simplificar nuestro Hamiltoniano si pasamos los operadores, tal y como hemos hecho en el anterior cap´ıtulo, al espacio Fourier como sigue1: ak=1 √N N−1 X j=0 eikjaj−→ aj=1 √N N−1 X k=0 e−ikj akbk=1 √N N−1 X j=0 eikjbj−→ bj=1 √N N−1 X k=0 e−ikjbk(4.2) a† k=1 √N N−1 X j=0 e−ikj a† j−→ a† j=1 √N N−1 X k=0 eikja† kb† k=1 √N N−1 X j=0 e−ikj b† j−→ b† j=1 √N N−1 X k=0 eikjb† k(4.3) Sustituyendo estas ´ultimas expresiones en la ecuaci´on (4.1) y haciendo uso de las relaciones de conmutaci´on ya demostradas [ai, a† j]=[bi, b† j] = δij [ai, b† j] = 0, obtenemos la expresi´on: H= N−1 X k=0 ωa† kak+b† kbk(ω+ 2Jcos k) + J0(a† kbk(1 + e−ik) + b† kak(1 + eik)) Para diagonalizar la matriz podemos observar que si definimos los vectores: ~ck=ak bk~c† k=a† kb† k Podemos expresarlo como: H= N−1 X k=0 ~c† kω0 0ω+ 2Jcos k~ck+~c† k0J0(1 + e−ik) J0(1 + eik) 0 ~ck O lo que es lo mismo: H= N−1 X k=0 ~c† kω J0(1 + e−ik) J0(1 + eik)ω+ 2Jcos k~ck(4.4) Es decir es una matriz diagonal por cajas, si diagonalizamos cada caja:  ω− J0(1 + e−ik) J0(1 + eik) (ω+ 2Jcos k)−= 0 (ω−)(ω+ 2Jcos k−)−2J02(1 + cos k)=0 2−(2ω+ 2Jcos k) + 2 cos k(Jω −J02)+(ω2−2J02) = 0 Resolviendo la ecuaci´on de segundo grado obtenemos finalmente las bandas de energ´ıa buscadas: =Jcos k+ω±qJ2cos2k+ 2J02(cos k+ 1) (4.5) Analizaremos la existencia de bandas planas en las siguientes secciones (4.3 y4.4 ) en las que consideraremos un caso m´as general en que las frecuencias de oscilaci´on son distintas para los osciladores de la fila de arriba en relaci´on con los de abajo en la topolog´ıa de diente de sierra. Las bandas planas son aquellas en las que como su nombre indica son niveles de energ´ıa “planos” en los que por tanto se caracterizan por una velocidad de grupo nula ( d dk = 0). Adem´as en el caso de la cadena de infinitos osciladores estas bandas se caracterizan por estar infinitamente degeneradas. Para ver los c´alculos realizados en esta secci´on en m´as detalle cons´ultese el Ap´endice B al final del escrito. 1Haciendo iguales consideraciones que cap´ıtulo anterior en lo que respecta al par´ametro de red con una primera zona de Brillouin k∈[−π, π] (v´ease la p´agina 8). 12 4.2. Diente de sierra, caso de frecuencias distintas. Por generalizar el modelo consideraremos ahora que las frecuencias de los osciladores de la fila de arriba2tienen una frecuencia ωay los osciladores de la fila de abajo tienen una frecuencia ωb. Con estos retoques el Hamiltoniano queda: H= N−1 X j=0 ωaa† jaj+ωbb† jbj+J N−1 X j=0 b† jbj+1 +b† j+1bj+J0 N−1 X j=0 a† jbj+a† jbj+1 +b† jaj+b† j+1aj(4.6) Igual que en el caso anterior buscaremos obtener las bandas de energ´ıa y como se trata de un modelo Tight Binding aplicaremos un m´etodo an´alogo a la secci´on anterior. Para ello empezaremos pasando los operadores de creaci´on y destrucci´on al espacio Fourier: ak=1 √N N−1 X j=0 eikjaj−→ aj=1 √N N−1 X k=0 e−ikjakbk=1 √N N−1 X j=0 eikjbj−→ bj=1 √N N−1 X k=0 e−ikjbk(4.7) a† k=1 √N N−1 X j=0 e−ikja† j−→ a† j=1 √N N−1 X k=0 eikja† kb† k=1 √N N−1 X j=0 e−ikjb† j−→ b† j=1 √N N−1 X k=0 eikjb† k(4.8) Sustituyendo estas ´ultimas expresiones en la ecuaci´on (4.1) y haciendo uso de las relaciones de conmutaci´on, obtenemos la expresi´on: H= N−1 X k=0 ωaa† kak+b† kbk(ωb+ 2Jcos k) + J0(a† kbk(1 + e−ik) + b† kak(1 + eik)) Para diagonalizar la matriz podemos observar que si definimos los vectores: ~ck=ak bk~c† k=a† kb† k Podemos expresarlo como: H= N−1 X k=0 ~c† kωa0 0ωb+ 2Jcos k~ck+~c† k0J0(1 + e−ik) J0(1 + eik) 0 ~ck O lo que es lo mismo: H= N−1 X k=0 ~c† kωaJ0(1 + e−ik) J0(1 + eik)ωb+ 2Jcos k~ck(4.9) Es decir es una matriz diagonal por cajas, si diagonalizamos cada caja:  ωa− J0(1 + e−ik) J0(1 + eik) (ωb+ 2Jcos k)−= 0 (ωa−)(ωb+ 2Jcos k−)−2J02(1 + cos k) = 0 2−((ωa+ωb)+2Jcos k) + 2 cos k(Jωa−J02)+(ωaωb−2J02) = 0 Resolviendo la ecuaci´on de segundo grado: =Jcos k+ωa+ωb 2±sJcos k+ωa+ωb 22 −2 cos k(Jωa−J02)−(ωaωb−2J02) (4.10) 2V´ease figura 4.1, p´agina 11. 13 Podemos obtener una ecuaci´on m´as sencilla si lo podemos en funci´on de la frecuencia media Ω: Ω = ωa+ωb 2→ωb= 2Ω −ωa Obtenemos as´ı finalmente una versi´on simplificada de las bandas de energ´ıa: =Jcos k+ Ω ±q(Jcos k+ (Ω −ωa))2+ 2J02(1 + cos k)) (4.11) Para ver los c´alculos realizados en m´as detalle cons´ultese el Ap´endice C al final del escrito. Podemos a modo de prueba representar las bandas para Ω = 2, J0= 1, J= 2 y ωa= 3 y ver as´ı su aspecto: Figura 4.2: Bandas de energ´ıa Ω = 2, J0= 1, J= 2 y ωa= 3 con +→azul y −→rojo 14 4.3. C´alculo de las bandas planas para dos frecuencias distintas Como necesitamos comprobar que el c´alculo num´erico que desarrollaremos en el cap´ıtulo siguiente obtendremos las condiciones que hay que imponer a la ecuaci´on anterior (4.11) para que aparezcan bandas planas y as´ı precisar de una manera m´as de comprobar que nuestro c´alculo num´erico es correcto. Adem´as desarrollaremos por casos las distintas condiciones que tienen que cumplirse para que aparezcan bandas planas y ver la dependencia de la energ´ıa con los distintos par´ametros del sistema. Por tanto partiendo de la ´ultima expresi´on: ±=Jcos k+ ∆Ω + ωa±q(Jcos k+ ∆Ω)2+ 2J02(1 + cos k) (4.12) En donde hemos expresado la ecuaci´on anterior (4.11) en funci´on de ∆Ω = Ω −ωa, para buscar las bandas planas imponemos que la velocidad de grupo es nula es decir: d± dk =−Jsin k∓J02sin k+Jsin k(Jcos k+ ∆Ω) q(Jcos k+ ∆Ω)2+ 2J02(1 + cos k) = 0 (4.13) Procedemos a despejar J0en funci´on de J: J2[(Jcos k+ ∆Ω)2+ 2J02(1 + cos k)] = J2(Jcos k+ ∆Ω)2+J04+ 2JJ02(Jcos k+ ∆Ω) 2J02J2(1 + cos k) = J04+ 2J2J02(1 + cos k)−2J2J02+ 2JJ02∆Ω J04−2J02J2+ 2J02J∆Ω = 0 Obteniendo finalmente para casos J06= 0, la relaci´on: J0=±p2J(J−∆Ω) (4.14) Notemos una cosa: J04−2J02J2+ 2J02J∆Ω = 0 J02(J02−2J2+ 2J∆Ω) = 0 J0 J2 = 21−∆Ω J J0=Js21−∆Ω J(4.15) Luego si J > 0 entonces J0>0 y viceversa, luego esta es en realidad la ecuaci´on correcta. De hecho si se considera que el signo de J0es arbitrario y no depende de Jentonces cuando se procede a usar la f´ormula 6 del art´ıculo [6] uno encuentra problemas al calcular as´ı los estados localizados. Podemos obtener el mismo resultado (Ecuaci´on (4.14) y (4.15)) con otros m´etodos alternativos de c´alculo para ello v´ease el Ap´endice D. 15 Notemos que si las dos frecuencias son iguales (∆Ω = 0) entonces: J0=√2Jes la condici´on para que aparezca una banda plana. Se tiene que cumplir 2J(J−∆Ω) ⩾0 si y solo si: Si: J⩾0−→ J⩾∆Ω Si: J⩽0−→ J⩽∆Ω Sustituyendo (4.15) en (4.14) tenemos: ±=Jcos k+ Ω± | Jcos k+ 2J−∆Ω |(4.16) Veamos los distintos casos: Caso 1: Si Jcos k+ 2J−∆Ω >0, obtenemos como soluciones: a)+= 2Jcos k+ 2J+ωa b)−=−2J+ωb Si llamamos ξ≡2J−∆Ω. Entonces: Jcos k+ξ > 0, imponiendo para que se cumpla ∀ktenemos: Caso de J > 0:Si ξ > J −→ 2J−∆Ω > J −→ J > ∆Ω Caso de J < 0:Si ξ > −J−→ 2J−∆Ω >−J−→ J > ∆Ω 3 Caso 2: Si Jcos k+ 2J−∆Ω <0, obtenemos como soluciones: a)+=−2J+ωb b)−= 2Jcos k+ 2J+ωa Si llamamos ξ≡2J−∆Ω. Entonces: Jcos k+ξ < 0, imponiendo para que se cumpla ∀ktenemos: Caso de J > 0:Si ξ < J −→ 2J−∆Ω < J −→ J < ∆Ω Caso de J < 0:Si ξ < −J−→ 2J−∆Ω <−J−→ J < ∆Ω 3 16 Recordando que se ten´ıa que cumplir que: 2J(J−∆Ω) ⩾0, tenemos: Caso A: Se tiene que cumplir una de estas condiciones siguientes: J∈[0,+∞) si ∆Ω <0 J∈[∆Ω,+∞) si ∆Ω >0 Tenemos dos casos posibles seg´un lo dicho anteriormente: Caso 1 de J > 0: Entonces J > ∆Ω −→ (+= 2Jcos k+ 2J+ωa −=−2J+ωb Por ejemplo si tomamos J= 5, ωa= 6 y un valor de ∆Ω que nos garantice que ωb>0 y que se cumple que J < ∆Ω por ejemplo ∆Ω = −2 con estos datos tenemos: (+= 10 cos k+ 16 −=−8 Si lo representamos en la primera zona de Brillouin pero usando directamente la ecuaci´on (4.11) directamente con los datos y condiciones impuestos obtenemos: Figura 4.3: Bandas de energ´ıa J= 5,∆Ω = −2y ωa= 6 con +→azul y −→rojo Caso 2 de J > 0:Entonces J < ∆Ω −→ Caso imposible 17 Caso B: Se tiene que cumplir una de estas condiciones siguientes: J∈(−∞,∆Ω] si ∆Ω <0 J∈(−∞,0] si ∆Ω >0 Tenemos dos casos posibles seg´un lo dicho anteriormente: Caso 1 de J < 0:Entonces J > ∆Ω/3−→ Caso imposible Caso 2 de J < 0:Entonces J < ∆Ω/3−→ (+=−2J+ωb −= 2Jcos k+ 2J+ωa Por ejemplo si tomamos J=−5, ωa= 6 y un valor de ∆Ω que nos garantice que ωb>0 y que se cumple que J < ∆Ω/3 por ejemplo ∆Ω = −2 con estos datos tenemos: (+= 12 −=−10 cos k−4 Si lo representamos en la primera zona de Brillouin pero usando directamente la ecuaci´on (4.11) directamente con los datos y condiciones impuestos obtenemos: Figura 4.4: Bandas de energ´ıa J=−5,∆Ω = −2y ωa= 6 con +→azul y −→rojo 18 Debemos darnos cuenta que si no se cumplen las condiciones para Jde los casos Caso A yCaso B entonces no existir´a un J0ya que 2J(J−∆Ω) <0 y no existir´a la ra´ız, por lo tanto la introducci´on de dos frecuencias distintas implica que han de cumplirse unas condiciones que restrinjan Jpara que puedan aparecer ondas planas. Podemos hacer una representaci´on conjunta del Caso A y el Caso B en tres dimensiones representando las ecuaciones (4.12) imponiendo la condici´on (4.15) para valores fijos de ∆Ω = −2 y ωa= 6 en funci´on de k∈[−π, π] y J∈[−10,10], as´ı podemos ver la evoluci´on de las bandas. Sin embargo hay un intervalo de la gr´afica que no es v´alido cuando J∈[−2,0] ya que en ese intervalo 2J(J−∆Ω) <0 y J0ser´ıa complejo y no existir´ıa esa soluci´on. Por tanto sabiendo esto la representaci´on ser´ıa: Figura 4.5: Bandas de energ´ıa J∈[−10,10],∆Ω = −2y ωa= 6 con +→rojo y −→azul 19 As´ı, si representamos todo esto: Figura 5.3: Bandas de energ´ıa obtenidas num´ericamente J= 5,∆Ω = −2y ωa= 6 con L´ınea Verde −→ Nivel de energ´ıa m´aximo de la banda de energ´ıa superior, L´ınea Rosa −→ Nivel de energ´ıa m´ınimo de la banda de energ´ıa superior, L´ıneas Azul y Rojo −→ niveles de energ´ıa intermedios que aparecen debido a condiciones de frontera libres yL´ınea Negra −→ Nivel de energ´ıa de la banda plana. Podemos ahora representar el “gap” (“hueco”) de energ´ıa entre las dos bandas. Definido como la diferencia de energ´ıas entre el nivel m´as inferior de la banda superior y el valor de la banda plana, observando nuevamente que tiende al valor esperado que en este caso es 14: Figura 5.4: Gap de energ´ıa J= 5,∆Ω = −2y ωa= 6 Por tanto podemos decir que los c´alculos num´ericos dependen del tama˜no de la red as´ı a mayor n´umero de sitios m´as se asemeja el c´alculo te´orico al num´erico tendiendo al valor esperado para el tama˜no del gap, es decir observamos los efectos del numero finito de sitios. 26 5.2. M´etodo num´erico 1 Empecemos con el primer m´etodo para obtener la matriz del Hamiltoniano para un n´umero dado de excitaciones4 y de sitios. Usaremos el programa Mathematica para programar el m´etodo en cuesti´on. En primer lugar pensamos este m´etodo que hemos llamado 1 pero podr´ıamos llamarlo trivial en el sentido de que no presenta ninguna complicaci´on pero sin embargo este es muy ineficiente. Para este m´etodo debemos recordar el aspecto de las matrices de los operadores de creaci´on y destrucci´on que tienen una dimensi´on “n+1” donde “n” es el n´umero de excitaciones: c=           0√1 0 0 . . . . . . . . . 0 0 √2 0 . . . . . . . . . 0 0 0 √3. . . . . . . . . . . .. . .. . .. . .... 0 0 0 0 0 √n . . . . . .. . .. . .. . .. . .. . .           c†=             0 0 0 0 . . . . . . . . . √1 0 0 0 . . . . . . . . . 0√2 0 0 . . . . . . . . . 0 0 √3 0 . . . . . . . . . . . .. . .. . .. . .... 0 0 0 0 √n+ 1 0 . . . . . .. . .. . .. . .. . .. . .             (5.5) Lo que haremos es lo siguiente, definiremos cada operador de creaci´on de cada sitio (cada uno con su espacio) y luego obtendremos su equivalente en el espacio completo de todos los sitios multiplicando tensorialmente cada operador. Por ejemplo si estamos en un caso con 5 nodos de la red de diente de sierra tenemos que definir un operador de creaci´on y otro de destrucci´on para cada sitio en el espacio completo, es decir, el producto tensorial de los cinco espacios, por ejemplo el operador c† 3del espacio completo ser´ıa5:c† 3=I(1) ×I(2) ×c† 3(3) ×I(4) ×I(5). Una vez tenemos todos los operadores definidos los multiplicamos y sumamos de acuerdo con la expresi´on del Hamiltoniano (5.3). Podr´ıa parecer que ya tenemos la matriz del Hamiltoniano, sin embargo el producto tensorial de las bases de cada uno de los espacios no es la base del espacio de los estados de “n” excitaciones que es la que estamos buscando. De hecho la base del espacio6del sitio jen el caso de “n” excitaciones es Aj n≡ {|0i,|1i,|2i, . . . , |n−1i,|ni} y si tenemos Nsitios el espacio del operador creaci´on/destrucci´on del espacio completo es el producto tensorial de espacios es decir B≡A1 n×A2 n×···×AN−1 n×AN n | {z } Nveces y este no es el espacio de los estados de “n” excitaciones ya que por ejemplo volviendo al ejemplo de 5 sitios y 3 excitaciones la base Bcontiene |0,1,2,0,0ipero tambi´en contiene estados como |0,3,3,2,0ien el que la red tiene m´as de 3 excitaciones7. Por ello deberemos primero encontrar la base que proyecte el Hamiltoniano que hemos obtenido en el espacio de los estados de “n” excitaciones. Deberemos obtener la base de ese espacio. Para ello podemos darnos cuenta del siguiente hecho, los elementos de la base del espacio completo B≡ {|0,0, . . . , 0,0i, . . . , |n, n, . . . , n, ni} pueden ordenarse como si fuesen una sucesi´on de n´umeros de N cifras y en la base num´erica n+ 1 desde el n´umero 00 . . . 00 hasta el nn . . . nn y de estos seleccionar aquellos cuya suma de cifras sea “n” as´ı obtendremos los estados de “n” excitaciones. Notar que de esta manera cada autoestado tiene asociado un n´umero ´unico. As´ı por ejemplo el estado |0,1,1,1,0i que en el caso de tres excitaciones equivale al n´umero en base 4: 01110 o el n´umero decimal 44 que es exclusivo de ese estado, veremos que este hecho nos ser´a ´util para el segundo m´etodo num´erico m´as eficiente. Con la base obtenida se calcula la matriz de proyecci´on y multiplicando el Hamiltoniano por esta se obtiene la matriz del Hamiltoniano definitiva y de esta ya podemos obtener los autovalores y autovectores. Sin embargo tenemos un problema que la dimensi´on de la matriz crece muy deprisa de hecho si tenemos “n” excitaciones y “N” sitios crece como (n+ 1)Nque por ejemplo para el caso de 3 excitaciones y 7 sitios es (3 + 1)7= 16384 un tama˜no por lo menos respetable. Esto nos impide llegar a caso con muchos sitios lo que nos impide trabajar en sistemas con un gran n´umero de sitios no pudiendo observar efectos de tama˜no finito, as´ı como una gran lentitud en los c´alculos num´ericos. Para ver una explicaci´on m´as detallada del programa desarrollado en Mathematica cons´ultese el Ap´endice E. 4M´as de una excitaci´on. 5Con: I(i) matriz identidad del espacio i-´esimo. 6Es decir el espacio de las matrices (5.5) 7En concreto 8. 27 5.3. M´etodo num´erico 2 5.3.1. Construyendo la base: Como hemos dicho el problema est´a en la dimensi´on de la matriz del Hamiltoniano que crece demasiado deprisa para nuestros prop´ositos. Es por ello que precisamos de un m´etodo que reduzca dicho tama˜no y as´ı reduzca los tiempos de c´alculo. El m´etodo alternativo consistir´a en calcular la matriz del Hamiltoniano directamente en la base de los autovectores con “n” excitaciones. Para construir la matriz necesitamos construir la base, por ello hacemos uso de lo expuesto en la secci´on anterior. Como en el conjunto de todos los estados de la base completa Bse encuentran los estados de la base buscada (es decir es un subespacio de B) podemos construir la base completa que va desde el estado |0,0, . . . , 0,0ihasta el |n, n, . . . , n, niy de ella seleccionar aquellos estados que contengan solo “n” excitaciones. Para construir la base completa podemos recurrir a la propiedad de que esta puede ser ordenada como si cada estado fuese un n´umero de “N” cifras (si hay “N” sitios) en la base num´erica “n+1” (si hay n excitaciones). As´ı el estado |0,0, . . . , 0,1iequivale al n´umero 00 . . . 01 de manera que la base completa es la sucesi´on de n´umeros en base “n+1” de “N” cifras en orden ascendente desde el 00 . . . 00 hasta el nn . . . nn. Como Mathematica permite generar tanto n´umeros como vectores de “N” componentes que sean un n´umero en una base determinada generando esa sucesi´on de vectores cuyas componentes son como las cifras de los n´umeros de la sucesi´on podemos seleccionar aquellos vectores cuya suma de componentes sea “n” y tendremos as´ı la base construida. De esta manera la posici´on del n´umero en la sucesi´on inicial desde el 00 . . . 00 hasta el nn . . . nn tiene asignado un n´umero que se corresponde con su posici´on en la lista menos uno. As´ı por ejemplo para una excitaci´on base 2 el |0,0, . . . , 0,0i ≡ 00 . . . 00 equivale en decimal 0 y suposici´on en la sucesi´on es la 1 ´o |0,0, . . . , 0,1,0i ≡ 00 . . . 010 equivale al n´umero en decimal al 2 y su posici´on es la 3 en la sucesi´on, de manera que cada estado tiene asociado un n´umero propio y podemos “transformar” un estado en otro simplemente sumando la cantidad adecuada en cada caso. Esta propiedad la usaremos m´as adelante. Concluir´e esta subsecci´on con un ejemplo consideremos el caso de 5 sitios y una excitaci´on (base 2), primero construimos la base completa (sucesi´on inicial):             |0,0,0i |0,0,1i |0,1,0i |0,1,1i |1,0,0i |1,0,1i |1,1,0i |1,1,1i             Seguidamente seleccionamos los estados que tengan 1 excitaci´on y creamos la sucesi´on complementaria de los n´umero decimales a los cuales equivalen:   1 2 4 ↔  (0,0,1) (0,1,0) (1,0,0)  (5.6) As´ı tenemos construida la base buscada. 28 5.3.2. Construyendo la matriz: Una vez que ya tenemos la base buscada tanto como una sucesi´on de vectores como una sucesi´on de n´umeros tenemos ahora que construir la matriz. Para comenzar deberemos ver como act´uan los operadores de creaci´on y destrucci´on recordemos que: c|ni=√n|n−1i(5.7) c†|ni=√n+ 1 |n+ 1i(5.8) Es decir que por ejemplo si: c† 1|1,1,0i=√2|2,1,0i´o c2|1,1,0i=√1|1,0,0i. En el m´etodo num´erico modificaremos el estado | i (un estado cualquiera) y luego multiplicaremos por el coeficiente √(un coeficiente cualquiera, por ejemplo en (5.8) es √n+ 1). Podemos as´ı darnos cuenta que el operador c† 1en nuestro caso equivale a sumar el n´umero en base 3: 100 al n´umero en base 3 (estado de partida) 110 obteniendo 100 + 110 = 210 as´ı mismos el caso de c2equivale a 110 −010 = 100, resta puesto que es el operador destrucci´on. Ambos casos para simplificar los c´alculos y no tener que sumar o restar vectores podemos hacerlo en su equivalente decimal as´ı: 9 + 12 = 21 y 12 −3 = 9. Es decir los operadores equivalen a los n´umeros: c† 1≡9 y c2≡ −3. Por ello que los operadores c†ycequivalen a las listas de n´umeros que representan todos los operadores creaci´on y destrucci´on respectivamente, as´ı en nuestro ejemplo de 3 sitios y una excitaci´on (base 2) ser´ıan: c†≡  4 2 1 c≡ −4 −2 −1 (5.9) As´ı para el caso de c†, tenemos que c† 1≡4c† 2≡2 y c† 3≡1. Es decir los operadores de creaci´on y destrucci´on est´an asociados a las potencias ±(n+1)jdonde “n+1” es la base num´erica de los n´umeros de las sucesiones que coincide con el n´umero de excitaciones menos uno y jequivale a la posici´on de la cifra (o componente si es un vector) que indica la posici´on del sitio en el arreglo usado y el signo viene a ser positivo si crea una excitaci´on o negativo si la destruye. Para poder construir la matriz debemos recordar la forma del Hamiltoniano (luego veremos el caso no lineal): H= N X j=1 ωac† 2jc2j+ωbc† 2j−1c2j−1+J N X j=1 c† 2j−1c2j+1 +c† 2j+1c2j−1+J0 N X j=1 c† 2jc2j−1+c† 2jc2j+1 +c† 2j−1c2j+c† 2j+1c2j Para obtener la matriz del Hamiltoniano en la base del subespacio de “n” excitaciones lo que haremos ser´a dividirlo en la suma de operadores m´as sencillos, consideraremos cada pareja de operadores de creaci´on y destrucci´on como un ´unico operador: c† 2jc2j, c† 2j−1c2j−1, c† 2j−1c2j+1, c† 2j+1c2j−1, c† 2jc2j−1, c† 2jc2j+1, c† 2j−1c2jyc† 2j+1c2j. Para obtener la matriz del Hamiltoniano obtendremos todas las im´agenes de los elementos de la base, para ello haremos actuar los operadores anteriores sobre cada elemento de la base sumando a cada componente del vector (5.6) las correspondientes componentes de los vectores c†yc(Eq. (5.9)). Como estos operadores conservan el n´umero de excitaciones la suma de estos tres n´umeros (c†,cy (5.6)) ser´a un n´umero que este en la lista dada por (5.6) y por tanto un vector de la base. Pongamos un ejemplo el operador c† 1c3, para el caso de 3 sitios y 1 excitaci´on, que actuar´a sobre los elementos de la base en (5.6) {|0,0,1i,|0,1,0i,|1,0,0i} que equivalen a los n´umeros (los estado son como n´umero binarios) (v´ease (5.6)) 1,2,4. Los operadores c† 1yc3equivalen a sumar 4 y sumar -1 (haciendo uso de (5.9)). As´ı c† 1c3equivale a sumar 3. De manera que las im´agenes de los vectores de la base 1,2,4 son 4,5,7. Notar que 4 si que es un vector de la base pero no lo son 5 y 7. Notemos que las im´agenes de 2 y 4 ´o |0,1,0iy|1,0,0ison cero ya que no tienen excitaciones en el tercer sitio. Usaremos este m´etodo con los seis ´ultimos operadores 8ya que los dos primeros (c† 2jc2jyc† 2j−1c2j−1) son matrices diagonales como veremos. Una vez hallamos obtenido las im´agenes de todos los vectores de la base sobre los seis tipos de operadores para cada valor de j, para cada sitio, ya podr´ıamos construir las matrices. Comparando las im´agenes de los vectores de la base con los vectores de la base. As´ı en el ejemplo anterior compararemos 1,2,4 con 4,5,7. De esta manera podemos construir la matriz del operador c† 1c3. Comparando la primera imagen 4 con la base colocaremos en la primera columna de la matriz de c† 1c31 si coincide con la imagen y 0 si no por tanto la primera columna es (0 0 1)T (en vertical). El resto de im´agenes 5 y 7 no coinciden con ning´un vector de la base original y las columnas de la matriz ser´an (0 0 0)TDe esta manera quedan solucionados los problemas del tipo que c3|0,1,0i= 0 y con nuestro m´etodo: 010 −001 = 001 ↔2−1 = 1. 8c† 2j−1c2j+1, c† 2j+1c2j−1, c† 2jc2j−1, c† 2jc2j+1, c† 2j−1c2jyc† 2j+1c2j 29 Por tanto la matriz del ejemplo c† 1c3queda: c† 1c3≡  000 000 100 (5.10) Solo nos falta calcular los factores √×√, que aparecen al hacer actuar un operador destrucci´on ciy seguidamente operador creaci´on c† j, usando la sucesi´on de los vectores iniciales (5.6) antes de obtener sus im´agenes y sabiendo como act´uan los operadores (5.8) y (5.9). De manera que, volviendo al ejemplo, usando ahora la lista de vectores (que no de los n´umeros 1, 2 y 4):   (0,0,1) (0,1,0) (1,0,0)  (5.11) Usando (5.8) y (5.9) calculamos los factores √×√del ejemplo anterior usando (5.11), al aplicar c† 1c3sobre cada elemento de la base, obteniendo:   √0+1√1 √0+1√0 √1+1√0 =  1 0 0 (5.12) Multiplicamos cada columna de la matriz (5.10) por su coeficiente √×√, en este caso queda igual pero no siempre ocurre esto. As´ı, la primera columna la multiplicar´ıamos por el primer elemento de la lista anterior 1 y las otras dos por 0. Si se diera el caso de que, al cambiar los operadores de creaci´on y destrucci´on por sumas o restas de n´umeros dieran una imagen que no es correcta al aplicarlo sobre un estado de la base. Como cuando el operador destrucci´on act´ua sobre un sitio donde no hay excitaciones, entonces el factor √×√es cero por (5.7). As´ı construimos las matrices a partir de los estados de la base compar´andolos con sus im´agenes (notar que compararemos n´umeros decimales entre s´ı ya que hemos reducido tanto la base como las im´agenes de la misma a dos listas de n´umeros decimales) y multiplicaremos cada elemento de la matriz reci´en construida por su factor correspondiente √×√. As´ı repitiendo el proceso para todo jobtenemos todos los operadores que como hemos dicho son las parejas de operadores creaci´on y destrucci´on. Solo quedan las matrices de los operadores c† 2jc2jy c† 2j−1c2j−1, que son dos matrices diagonales. En la que en la primera los elementos impares de la diagonal son nulos y viceversa en la segunda con sus correspondientes factores √×√. El Hamiltoniano lo construimos finalmente multiplicando y sumando todas las matrices (las que representan a las parejas: c† 2jc2j, c† 2j−1c2j−1, c† 2j−1c2j+1, c† 2j+1c2j−1, c† 2jc2j−1, c† 2jc2j+1, c† 2j−1c2jyc† 2j+1c2j) y a continuaci´on podemos obtener con la matriz ya construida los autovalores y autovectores, matriz que no tenemos que proyectar porque la hemos construido en la base que nos interesaba. Obtendremos ahora la dimensi´on de la matriz para compararla con el m´etodo anterior y ver si realmente hemos simplificado el c´alculo. Para ello necesitamos calcular el n´umero de elementos de la base dado un “N” y un “n”. El problema se puede simplificar a el n´umero de combinaciones de meter “n” bolas (excitaciones) en “N” cajas (“sitios”), o lo que es lo mismo ordenar (“N-1”) “separadores” y n “bolas” es decir (N−1+n)!, pero como podemos intercambiar dos bolas o separadores entre s´ı sin cambiar de estado tenemos entonces que quitar las permutaciones de estos obteniendo: (N−1 + n)! n!(N−1)! =N−1 + n n Que si lo comparamos con (n+ 1)Nse puede demostrar que N−1 + n n<(n+ 1)Nde hecho para n= 5 y N= 7 tenemos que 7−1+5 5= 462 (5+1)7= 279936. Esto solo es el tama˜no del lado de la matriz cuadrada, el total de n´umeros de cada matriz corresponde al cuadrado de esas cantidades. Podemos ver en detalle explicado el programa hecho para Mathematica y para ello cons´ultese el Ap´endice F. 30 5.4. Comprobaci´on de los m´etodos num´ericos Un primer chequeo que podemos hacer para comprobar que el nuevo programa calcula correctamente la matriz del Hamiltoniano comparamos los valores de los autovalores en ambos casos. Usaremos como par´ametros ωa= 2, J= 5, ∆Ω = −2 y 5 sitios compraremos para distinto n´umero de excitaciones n= 1,2,3. Con esto en las tablas 5.2, 5.3 y 5.4 se puede notar, que la diferencia de los valores obtenidos es del orden del error num´erico: Cuadro 5.2: Caso de n= 1 excitaciones. M´etodo 2 M´etodo 1 DIFERENCIA ( %) 18,00 18,00 0,00E+00 8,60 8,60 0,00E+00 -8,00 -8,00 0,00E+00 -8,60 -8,60 1,24E-13 -12,00 -12,00 0,00E+00 Cuadro 5.3: Caso de n= 2 excitaciones. M´etodo 2 M´etodo 1 DIFERENCIA ( %) 36,00 36,00 2,76E-13 26,60 26,60 0,00E+00 17,20 17,20 0,00E+00 10,00 10,00 0,00E+00 9,40 9,40 9,45E-14 6,00 6,00 0,00E+00 0,60 0,60 1,99E-12 7,85E-15 -2,83E-16 103,61 -3,40 -3,40 3,01E-13 -16,00 -16,00 0,00E+00 -16,60 -16,60 0,00E+00 -17,20 -17,20 0,00E+00 -20,00 -20,00 0,00E+00 -20,60 -20,60 -4,83E-13 -24,00 -24,00 4,14E-13 31 Cuadro 5.4: Caso de n= 3 excitaciones. M´etodo 2 M´etodo 1 DIFERENCIA ( %) 54,00 54,00 -3,68E-13 44,60 44,60 0,00E+00 35,20 35,20 -3,03E-13 28,00 28,00 0,00E+00 27,40 27,40 3,76E-13 25,81 25,81 -3,99E-13 24,00 24,00 0,00E+00 18,60 18,60 0,00E+00 18,00 18,00 5,53E-13 14,60 14,60 0,00E+00 9,20 9,20 2,32E-13 8,60 8,60 -3,51E-13 5,20 5,20 0,00E+00 2,00 2,00 5,11E-13 1,40 1,40 0,00E+00 0,80 0,80 2,51E-13 -2,00 -2,00 0,00E+00 -2,60 -2,60 -3,92E-13 -6,00 -6,00 0,00E+00 -7,40 -7,40 -2,76E-13 -8,00 -8,00 0,00E+00 -8,60 -8,60 -1,24E-13 -11,40 -11,40 0,00E+00 -12,00 -12,00 -8,29E-13 -15,40 -15,40 0,00E+00 -24,00 -24,00 -4,14E-13 -24,60 -24,60 0,00E+00 -25,20 -25,20 0,00E+00 -25,81 -25,81 -3,99E-13 -28,00 -28,00 3,55E-13 -28,60 -28,60 0,00E+00 -29,20 -29,20 3,41E-13 -32,00 -32,00 3,11E-13 -32,60 -32,60 3,05E-13 -36,00 -36,00 0,00E+00 32 Otro m´etodo de comprobaci´on que aplicaremos es calcular con el M´etodo 2, los autovalores para el caso de “single excitation” (una sola excitaci´on) para los n´umero de sitios desde 5 hasta 17 y comprobar la dependencia del gap con el tama˜no de la red de osciladores dispuestos en forma de diente de sierra. Es decir que conforme se va aumentando el n´umero de sitios el tama˜no del gap tiende al valor te´orico, para un n´umero muy grande de sitios, asint´oticamente. As´ı tenemos que los autovalores son: Cuadro 5.5: Autovalores calculados a trav´es del M´etodo 2 para distinto n´umero de sitios y una excitaci´on. Para valores de ωa= 6, J= 5 y ∆Ω = −2. Aplicando las condiciones para J0para que exista una banda plana. Sitios Autovalores 5 -8,00 -4,60 -4,00 12,60 22,00 7 -8,00 -4,39 -4,24 9,84 17,24 23,56 9 -8,00 -4,34 -4,29 8,44 13,94 19,90 24,36 11 -8,00 -4,32 -4,31 7,66 11,77 16,84 21,54 24,82 13 -8,00 -4,32 -4,31 7,20 10,34 14,52 18,86 22,60 25,11 15 -8,00 -4,32 -4,32 6,90 9,35 12,80 16,63 20,31 23,33 25,31 17 -8,00 -4,32 -4,32 6,70 8,66 11,51 14,84 18,25 21,37 23,85 25,45 Si se obtienen las diferencias con los ya calculados en la primera secci´on de este cap´ıtulo que se calcularon construyendo la matriz directamente por que se sab´ıa cual era su forma. No se obtienen diferencias superiores a 10−13 lo cual verifica el correcto funcionamiento del m´etodo presentando exactamente la misma dependencia del gap con el tama˜no de la red. Este caso recordemos que corresponde al Caso A de la p´agina 17, lo que nos sirve de conexi´on entre el c´alculo anal´ıtico y num´erico. 33 5.5. Resultados 5.5.1. Caso lineal Entre los resultados que podemos comprobar primero para el caso no lineal se encuentra la f´ormula citada en la referencia [7] y [6] en la que nos dan una expresi´on anal´ıtica para obtener los autoestados localizados, es decir para el caso de una ´unica excitaci´on y en el caso de que las filas superior e inferior de la topolog´ıa en diente de sierra con frecuencias iguales, encontramos que los autoestados de la banda plana son: |Γji=1 2(c† 2j+c† 2j+2 −√2c† 2j+1)|0i(5.13) Notemos que la f´ormula dada en bibliograf´ıa [6] no es la expresi´on de arriba. Es una errata que se subsana consultando la bibliograf´ıa que proporciona ese art´ıculo que es [8]. Para comprobar el resultado con el programa es sencillo, no hay m´as que construir el estado y obtener su imagen haciendo actuar el Hamiltoniano reci´en obtenido, obtener la norma del vector imagen, dividirla por la norma del vector y obtener as´ı el valor del autovector. En donde aplicando estas modificaciones al M´etodo 2 podemos comprobar que la expresi´on (5.10) funciona. As´ı por ejemplo para un caso concreto de valores ωa= 6, J=−5 y ∆Ω = −2. Notar que los autoestados localizados tienen la forma en la topolog´ıa de: Figura 5.5: Un autoestado localizado. De manera que para este caso particular de nueve sitios la degeneraci´on de los autoestados localizados es de 3. Veamos la degeneraci´on del autoestado localizado para el caso de una excitaci´on. Veremos en la siguiente secci´on para el caso de m´as de una excitaci´on. Si “m” es la degeneraci´on el n´umero de sitios “N” que es un n´umero impar como ya hemos mencionado antes sigue: N= 2m+3 empezando para el caso de 5 sitios que es el m´ınimo n´umero de sitios en el cual tiene sentido que puedan aparecer estos autoestados localizados tal y como se puede observar en el dibujo. Es decir la degeneraci´on es: m=N−3 2en donde si obtenemos los autovalores para el caso de una excitaci´on y varios sitios con el M´etodo 2 como en el cuadro 5.5 obtenemos que la degeneraci´on del autoestado localizado9coincide exactamente con la expresi´on mencionada. As´ı en funci´on del n´umero de sitios (en el caso de una excitaci´on) la degeneraci´on del autoestado localizado para el cuadro 5.5 obtuvimos con el M´etodo 2, resultado v´alido solo para una excitaci´on10: Sitios Degeneraci´on 5 1 7 2 9 3 11 4 13 5 15 6 17 7 Cuadro 5.6: Degeneraci´on de los autoestados localizados obtenidos a trav´es del M´etodo 2 en funci´on del n´umero de sitios y una excitaci´on. Para valores de ωa= 6, J=−5 y ∆Ω = −2. Aplicando las condiciones para J0para que exista una banda plana. (Caso de la tabla anterior 5.5) Para ver cual es la modificaci´on que hay que a˜nadir al c´odigo del M´etodo 2 para comprobar que (5.13) es un autoestado localizado, cons´ultese el Ap´endice G. 9El autoestado cuyo autovalor es el de la banda plana en el caso anterior era =−8. 10En la siguiente secci´on veremos el caso general. 34 Nosotros vamos a intentar obtener una formula m´as general del caso anterior para los autoestados localizados para el caso de dos frecuencias distintas de la fila superior e inferior de sitios del arreglo de diente de sierra para el caso de una sola excitaci´on. Para ello veamos el ejemplo de 5 sitos y una excitaci´on que sabemos la forma de la matriz. Si aplicamos uno de los autoestados (5.13) para el caso de frecuencias iguales obtenemos:       ωbJ0J0 0 J0ωaJ00 0 J J0ωbJ0J 0 0 J0ωaJ0 0 0 J J0ωb             0 1/2 −√2/2 1/2 0       =1 2       J0−√2J ωa−√2J0 J0−√2ωb+J0 −√2J0+ωa −√2J+J0       (5.14) Que claramente no es autovector en el caso de frecuencias distintas como cabr´ıa esperar. Definamos tres constantes A,ByC, y veamos que valor podemos darles para que el anterior estado sea autovector del Hamiltoniano:       ωbJ0J0 0 J0ωaJ00 0 J J0ωbJ0J 0 0 J0ωaJ0 0 0 J J0ωb             0 A/2 −B√2/2 C/2 0       =1 2       AJ0−B√2J Aωa−B√2J0 AJ0−B√2ωb+CJ0 −B√2J0+Cωa −B√2J+CJ0       (5.15) Con ello imponemos: J0A−B√2J= 0 (5.16) J0C−B√2J= 0 (5.17) Es decir se tiene que cumplir que A=C. Con ello podemos hallar la relaci´on entre AyB: J0A−B√2J= 0 −→ A B=√2J J0=rJ J−∆Ω (5.18) En donde hemos hecho uso de la condici´on de banda plana: J0=Js21−∆Ω J. Por ´ultimo establecemos la condici´on para que los elementos del vector imagen sean proporcionales a los del vector inicial y as´ı sea autoestado: Aωa−B√2J0 A=Bωb−J0A+B √2 B(5.19) Sustituimos la condici´on (5.18) para que dicha relaci´on garantice que es un autoestado si se cumple (5.19), obteniendo: −√2J2+ (√2(ωb−ωa) + J0)J+√2J0 √2J= 0 Obteniendo la relaci´on entre JyJ0para que el estado sea autovector: J0=√2J(J−2∆Ω) J+√2(5.20) Pero al mismo tiempo queremos que existan bandas planas, ya que estamos buscando los autoestados localizados. Por ello a de cumplirse simult´aneamente: J0=Js21−∆Ω J(5.21) Igualando ambas expresiones: Js21−∆Ω J=√2J(J−2∆Ω) J+√2 35 Llegados a este punto la respuesta para poder continuar parec´ıa clara deb´ıa de dar con un m´etodo que nos permitiese obtener la base de ese subespacio, que conservaba el numero de excitaciones. Para as´ı reducir el tama˜no de nuestras matrices as´ı como del tiempo de c´alculo. Una vez obtenido procedimos a comparar ambos m´etodos, as´ı como con los c´alculos anal´ıticos, para comprobar que efectivamente calculaban correctamente la matriz del Hamiltoniano. Adem´as comprobamos usando una expresi´on para los autoestados localizados que encontramos en la bibliograf´ıa. Finalmente usamos el m´etodo m´as eficiente para calcular los niveles de energ´ıa con el caso no lineal. Para ello no tuvimos m´as que hacer una peque˜na modificaci´on al m´etodo para a˜nadir el t´ermino no lienal. Obtenidos los niveles de energ´ıa en funci´on de la constante no lineal de acoplamiento pudimos observar el desdoblamiento de energ´ıa, as´ı como obtener la expresi´on anal´ıtica que daba la degeneraci´on de los autoestados localizados. Para ello fue clave tanto conocer la forma de los autoestados, que nos proporcion´o la bibliograf´ıa, como el desdoblamiento que observamos de los niveles de energ´ıa. Esta vez hab´ıamos obtenido la degeneraci´on de los autoestados localizados para cualquier n´umero de excitaciones. Hemos as´ı obtenido un m´etodo num´erico que nos ha permitido resolver un caso no tratable num´ericamente, el caso no lineal. Las principales conclusiones son: Primero el c´alculo anal´ıtico de las bandas para el caso m´as general de frecuencias distintas nos ha llevado al la conclusi´on que estas solo existen bajo determinadas condiciones. En segundo lugar en cuanto al c´alculo num´erico del caso lineal ten´ıamos: El tama˜no finito del sistema es importante, es decir conforme aumentamos el tama˜no del sistema se aprecia una mejor concordancia por ejemplo en el c´alculo del gap con respecto al c´alculo anal´ıtico. En tercer lugar hemos concluido que deb´ıamos de llegar a un m´etodo que redujera los tiempos de c´alculo para poder abordar as´ı casos de m´as excitaciones y sitios. En cuarto lugar en el caso no lineal vimos el desdoblamiento de los niveles de energ´ıa degenerados en una banda plana para el caso lineal. Ello y el conocimiento previo de la “topolog´ıa” de los autoestados localizados no permiti´o obtener la degeneraci´on de los niveles de energ´ıa de la banda plana en funci´on del n´umero de sitios del arreglo y de las excitaciones. 42 Ap´endice A C´alculo de bandas de cadena infinita de osciladores acoplados Consideraremos el caso de una cadena infinita de osciladores bos´onicos sometidos a un potencial peri´odico acoplados a trav´es de Jque controla este acoplamiento entre primeros vecinos: Figura A.1: Cadena de osciladores. Vamos a obtener la bandas de energ´ıa del Hamiltoniano: H= N−1 X j=0 [ωa† jaj−J(aja† j+1 +a† jaj+1)] (A.1) Para ello haremos uso de las ecuaciones: ak=1 √N N−1 X j=0 eikjaj−→ aj=1 √N N−1 X k=0 e−ikjak(A.2) a† k=1 √N N−1 X j=0 e−ikja† j−→ a† j=1 √N N−1 X k=0 eikja† k(A.3) Adem´as de estas expresiones deberemos de hacer uso de las relaciones de conmutaci´on: [ak, a† k0] = δkk0[ai, a† j] = δij [ak, ak0] = 0 [a† k, a† k0] = 0 (A.4) 43 Empecemos el desarrollo sustituyendo las ecuaciones (A.1) y (A.3) en el Hamiltoniano: H=1 N N−1 X j=0 ω N−1 X kk0a† kak0ei(k−k0)j−J N−1 X kk0aka† k0ei(k0(j+1)−kj)+a† k0akei(k0j−k(j+1))= =1 NN−1 X kk0 ωa† kak0 N−1 X j=0 ei(k−k0)j−J N−1 X kk0 aka† k0eik0 N−1 X j=0 ei(k0−k)j−J N−1 X kk0 a† k0ake−ik N−1 X j=0 ei(k0−k)j= Puede demostrarse que PN−1 j=0 ei(k−k0)j=δ(k−k0). Haciendo uso de esta expresi´on podemos continuar el c´alculo anterior: =1 NN−1 X kk0 ωa† kak0Nδ(k−k0)−J N−1 X kk0 aka† k0eik0Nδ(k0−k)−J N−1 X kk0 a† k0ake−ikNδ(k0−k)= = N−1 X k=0 ωa† kak−J(aka† keik +a† kake−ik)= N−1 X k=0 ωa† kak−2Ja† kakeik +e−ik 2−J N−1 X k=0 e−ik = = N−1 X k=0 ωa† kak−2Ja† kakcos(k)−JNδ(1) = N−1 X k=0 (ω−2Jcos(k))a† kak= N−1 X k=0 ωka† kak Obteniendo as´ı la relaci´on buscada: H= N−1 X k=0 ωka† kak(A.5) 44 Ap´endice B C´alculo de las bandas en topolog´ıa de tipo dientes de sierra para frecuencias iguales Nuestro objetivo es encontrar las bandas de energ´ıa para el siguiente Hamiltoniano: H= N−1 X j=0 ω(a† jaj+b† jbj)+J N−1 X j=0 b† jbj+1 +b† j+1bj+J0 N−1 X j=0 a† jbj+a† jbj+1 +b† jaj+b† j+1aj(B.1) Que representa un conjunto de osciladores de frecuencia ωacoplados a trav´es de las constantes J0yJen un modelo Tight Binding para una topolog´ıa tal y como se muestra en la figura: Figura B.1: Disposici´on en diente de sierra de los sitios. Para obtener las bandas de energ´ıa en primer lugar pasaremos los operadores de creaci´on y destrucci´on al espacio Fourier usando para ello: ak=1 √N N−1 X j=0 eikjaj−→ aj=1 √N N−1 X k=0 e−ikjakbk=1 √N N−1 X j=0 eikjbj−→ bj=1 √N N−1 X k=0 e−ikjbk(B.2) a† k=1 √N N−1 X j=0 e−ikja† j−→ a† j=1 √N N−1 X k=0 eikja† kb† k=1 √N N−1 X j=0 e−ikjb† j−→ b† j=1 √N N−1 X k=0 eikjb† k(B.3) 45 Con todo esto sustituyendo las expresiones (B.2) y (B.3) en el Hamiltoniano (B.1) obtenemos, reorganizando los sumatorios: H=1 N ω N−1 X j=0 N−1 X kk0eij(k−k0)(a† kak0+b† kbk0)+J N−1 X j=0 N−1 X kk0eij(k−k0)(e−ikb† kbk0+eikb† kbk0)+ +J0 N−1 X j=0 N−1 X kk0eij(k−k0)(a† kbk0+b† kak0) + eij(k−k0)(a† kbk0e−ik0+b† kak0eik)! Reorganizando los sumatorios, tenemos: H=1 N ω N−1 X kk0 (a† kak0+b† kbk0)N−1 X j=0 eij(k−k0)+J N−1 X kk0 (e−ikb† kbk0+eikb† kbk0)N−1 X j=0 eij(k−k0)+ +J0 N−1 X kk0 (a† kbk0+b† kak0)N−1 X j=0 eij(k−k0)+J0 N−1 X kk0 (a† kbk0e−ik0+b† kak0eik)N−1 X j=0 eij(k−k0)! Haciendo uso de la propiedad N−1 X j=0 eij(k−k0)=Nδ(k−k0), tenemos: H=1 N ω N−1 X kk0 (a† kak0+b† kbk0)Nδ(k−k0) + J N−1 X kk0 (e−ikb† kbk0+eikb† kbk0)Nδ(k−k0)+ +J0 N−1 X kk0 (a† kbk0+b† kak0)Nδ(k−k0) + J0 N−1 X kk0 (a† kbk0e−ik0+b† kak0eik)Nδ(k−k0)! H=ω N−1 X k=0 (a† kak+b† kbk) + J N−1 X k=0 (e−ikb† kbk+eikb† kbk) + J0 N−1 X k=0 (a† kbk+b† kak) + J0 N−1 X k=0 (a† kbke−ik +b† kakeik) En donde usamos las relaciones de conmutaci´on: [ai, a† j] = [bi, b† j] = δij [ai, b† j] = 0 Continuando con el c´alculo: H=ω N−1 X k=0 (a† kak+b† kbk) + J N−1 X k=0 b† kbk2 cos k+J0 N−1 X k=0 (a† kbk+b† kak) + J0 N−1 X k=0 (a† kbke−ik +b† kakeik) H= N−1 X k=0 ωa† kak+b† kbk(ω+ 2Jcos k) + J0(a† kbk(1 + e−ik) + b† kak(1 + eik)) 46 Si ahora definimos los siguientes vectores: ~ck=ak bk~c† k=a† kb† k(B.4) Sustituyendo: H= N−1 X k=0 ~c† kω0 0ω+ 2Jcos k~ck+~c† k0J0(1 + e−ik) J0(1 + eik) 0 ~ck Obteniendo finalmente una matriz que es diagonal por cajas: H= N−1 X k=0 ~c† kω J0(1 + e−ik) J0(1 + eik)ω+ 2Jcos k~ck(B.5) Si no vemos claro que sea diagonal por cajas podemos definir estos vectores: ~ h=               . . . ak−1 bk−1 ak bk ak+1 bk+1 . . .               (B.6) ~ h†=. . . a† k−1b† k−1a† kb† ka† k+1 b† k+1 . . . (B.7) De manera que nos queda finalmente esta matriz: H=~ h†              .... . .. . .. . .. . .. . .. . . . . . ω J0(1 + e−i(k−1)) 0 0 0 0 . . . . . . J0(1 + ei(k−1))ω+ 2Jcos (k−1) 0 0 0 0 . . . . . . 0 0 ω J0(1 + e−ik) 0 0 . . . . . . 0 0 J0(1 + eik)ω+ 2Jcos (k) 0 0 . . . . . . 0 0 0 0 ω J0(1 + e−i(k+1)). . . . . . 0 0 0 0 J0(1 + ei(k+1))ω+ 2Jcos (k+ 1) . . . . . .. . .. . .. . .. . .. . ....              ~ h (B.8) Con todo ello si diagonalizamos cada caja, obtenemos:  ω− J0(1 + e−ik) J0(1 + eik) (ω+ 2Jcos k)−= 0 (ω−)(ω+ 2Jcos k−)−2J02(1 + cos k) = 0 2−(2ω+ 2Jcos k) + 2 cos k(Jω −J02)+(ω2−2J02)=0 Resolviendo la ecuaci´on de segundo grado, obtenemos finalmente las bandas de energ´ıas buscadas: =Jcos k+ω±qJ2cos2k+ 2J02(cos k+ 1) (B.9) 47 48 Ap´endice C C´alculo de las bandas en topolog´ıa de tipo dientes de sierra para frecuencias distintas Nuestro objetivo es encontrar las bandas de energ´ıa para el siguiente Hamiltoniano: H= N−1 X j=0 ωaa† jaj+ωbb† jbj+J N−1 X j=0 b† jbj+1 +b† j+1bj+J0 N−1 X j=0 a† jbj+a† jbj+1 +b† jaj+b† j+1aj(C.1) Que representa un conjunto de osciladores de frecuencia ωa(fila superior de osciladores [V´ease la siguiente figura]) yωb(fila inferior de osciladores) acoplados a trav´es de las constantes J0yJen un modelo Tight Binding para una topolog´ıa tal y como se muestra en la figura: Figura C.1: Disposici´on en diente de sierra de los sitios. Para obtener las bandas de energ´ıa en primer lugar pasaremos los operadores de creaci´on y destrucci´on al espacio Fourier usando para ello: ak=1 √N N−1 X j=0 eikjaj−→ aj=1 √N N−1 X k=0 e−ikjakbk=1 √N N−1 X j=0 eikjbj−→ bj=1 √N N−1 X k=0 e−ikjbk(C.2) a† k=1 √N N−1 X j=0 e−ikja† j−→ a† j=1 √N N−1 X k=0 eikja† kb† k=1 √N N−1 X j=0 e−ikjb† j−→ b† j=1 √N N−1 X k=0 eikjb† k(C.3) Con todo esto sustituyendo las expresiones (C.2) y (C.3) en el Hamiltoniano (C.1) obtenemos, reorganizando los sumatorios: 49 H=1 N N−1 X j=0 N−1 X kk0eij(k−k0)(ωaa† kak0+ωbb† kbk0)+J N−1 X j=0 N−1 X kk0eij(k−k0)(e−ikb† kbk0+eikb† kbk0)+ +J0 N−1 X j=0 N−1 X kk0eij(k−k0)(a† kbk0+b† kak0) + eij(k−k0)(a† kbk0e−ik0+b† kak0eik)! Reorganizando los sumatorios, tenemos: H=1 N N−1 X kk0 (ωaa† kak0+ωbb† kbk0)N−1 X j=0 eij(k−k0)+J N−1 X kk0 (e−ikb† kbk0+eikb† kbk0)N−1 X j=0 eij(k−k0)+ +J0 N−1 X kk0 (a† kbk0+b† kak0)N−1 X j=0 eij(k−k0)+J0 N−1 X kk0 (a† kbk0e−ik0+b† kak0eik)N−1 X j=0 eij(k−k0)! Haciendo uso de la propiedad N−1 X j=0 eij(k−k0)=Nδ(k−k0), tenemos: H=1 N N−1 X kk0 (ωaa† kak0+ωbb† kbk0)Nδ(k−k0) + J N−1 X kk0 (e−ikb† kbk0+eikb† kbk0)Nδ(k−k0)+ +J0 N−1 X kk0 (a† kbk0+b† kak0)Nδ(k−k0) + J0 N−1 X kk0 (a† kbk0e−ik0+b† kak0eik)Nδ(k−k0)! H= N−1 X k=0 (ωaa† kak+b† kbk) + J N−1 X k=0 (e−ikb† kbk+eikb† kbk) + J0 N−1 X k=0 (a† kbk+b† kak) + J0 N−1 X k=0 (a† kbke−ik +b† kakeik) En donde usamos las relaciones de conmutaci´on: [ai, a† j]=[bi, b† j] = δij [ai, b† j]=0 Continuando con el c´alculo: H= N−1 X k=0 (ωaa† kak+ωbb† kbk) + J N−1 X k=0 b† kbk2 cos k+J0 N−1 X k=0 (a† kbk+b† kak) + J0 N−1 X k=0 (a† kbke−ik +b† kakeik) H= N−1 X k=0 ωaa† kak+b† kbk(ωb+ 2Jcos k) + J0(a† kbk(1 + e−ik) + b† kak(1 + eik)) Si ahora definimos los siguientes vectores: ~ck=ak bk~c† k=a† kb† k(C.4) 50 Sustituyendo las anteriores expresiones: H= N−1 X k=0 ~c† kωa0 0ωb+ 2Jcos k~ck+~c† k0J0(1 + e−ik) J0(1 + eik) 0 ~ck Obteniendo finalmente una matriz que es diagonal por cajas: H= N−1 X k=0 ~c† kωaJ0(1 + e−ik) J0(1 + eik)ωb+ 2Jcos k~ck(C.5) Si no vemos claro que sea diagonal por cajas podemos definir estos vectores: ~ h=               . . . ak−1 bk−1 ak bk ak+1 bk+1 . . .               (C.6) ~ h†=. . . a† k−1b† k−1a† kb† ka† k+1 b† k+1 . . . (C.7) De manera que nos queda finalmente esta matriz: H=~ h†              .... . .. . .. . .. . .. . .. . . . . . ωaJ0(1 + e−i(k−1)) 0 0 0 0 . . . . . . J0(1 + ei(k−1))ωb+ 2Jcos (k−1) 0 0 0 0 . . . . . . 0 0 ωaJ0(1 + e−ik) 0 0 . . . . . . 0 0 J0(1 + eik)ωb+ 2Jcos (k) 0 0 . . . . . . 0 0 0 0 ωaJ0(1 + e−i(k+1)). . . . . . 0 0 0 0 J0(1 + ei(k+1))ωb+ 2Jcos (k+ 1) . . . . . .. . .. . .. . .. . .. . ....              ~ h (C.8) Con todo ello si diagonalizamos cada caja, obtenemos:  ωa− J0(1 + e−ik) J0(1 + eik) (ωb+ 2Jcos k)−= 0 (ωa−)(ωb+ 2Jcos k−)−2J02(1 + cos k) = 0 2−((ωa+ωb)+2Jcos k) + 2 cos k(Jωa−J02)+(ωaωb−2J02) = 0 Resolviendo la ecuaci´on de segundo grado, obtenemos finalmente las bandas de energ´ıas buscadas: =Jcos k+ωa+ωb 2±sJcos k+ωa+ωb 22 −2 cos k(Jωa−J02)−(ωaωb−2J02) (C.9) Podemos obtener una ecuaci´on m´as sencilla si lo podemos en funci´on de la frecuencia media Ω: Ω = ωa+ωb 2→ωb= 2Ω −ωa =Jcos k+ Ω ±q(Jcos k+ (Ω −ωa))2+ 2J02(1 + cos k)) (C.10) 51 E.2. Resumen del c´odigo El resumen de todo el c´odigo usado es: sitios = 5; dim = 4; n= dim −1; Id = IdentityMatrix[dim]; a= Table[0,{i, 1,dim},{j, 1,dim}]; For[i= 1, i < dim, i++, a[[i, i + 1]] = Sqrt[i]; ] adag = Transpose[a]; a1 = SparseArray[KroneckerProduct[a, Id,Id,Id,Id]]; a1d = SparseArray[KroneckerProduct[adag,Id,Id,Id,Id]]; a2 = SparseArray[KroneckerProduct[Id, a, Id,Id,Id]]; a2d = SparseArray[KroneckerProduct[Id,adag,Id,Id,Id]]; a3 = SparseArray[KroneckerProduct[Id,Id, a, Id,Id]]; a3d = SparseArray[KroneckerProduct[Id,Id,adag,Id,Id]]; a4 = SparseArray[KroneckerProduct[Id,Id,Id, a, Id]]; a4d = SparseArray[KroneckerProduct[Id,Id,Id,adag,Id]]; a5 = SparseArray[KroneckerProduct[Id,Id,Id,Id, a]]; a5d = SparseArray[KroneckerProduct[Id,Id,Id,Id,adag]]; Ω = −2.; (*Ω ≡∆Ω*) ω1=6.; ω2=2.∗Ω + ω1; J1 = −5; J2 = p2∗J1 ∗(J1 −Ω); (*El Hamiltoniano*) H= SparseArray[ω1∗(a2d.a2 + a4d.a4) + ω2∗(a1d.a1 + a3d.a3 + a5d.a5) + J1 ∗(a1d.a3 + a3d.a1 + a3d.a5 + a5d.a3) + J2 ∗(a2d.a1 + a4d.a3 + a2d.a3 + a4d.a5 + a1d.a2 + a3d.a4 + a3d.a2 + a5d.a4)]; Basis = Table IntegerDigits[m, dim,sitios],m, 0,(dim)sitios −1; Res = Table If[Sum[Basis[[j]][[i]],{i, 1,sitios}] == n, j, 0],j, 1,(dim)sitios; Res2 = Select[Res,#6= 0&]; Proyec = Table KroneckerDelta[m, i]∗KroneckerDelta[m, j]∗KroneckerDelta[i, j],{m, Res2},i, 1,(dim)sitios, j, 1,(dim)sitios; Proyec2 = Sum[Proyec[[m]][[]][[]],{m, 1,Length[Res2]}]; HProyec = Proyec2.H; Eigenvalues[HProyec] Eigenvectors[HProyec] 58 Ap´endice F M´etodo num´erico 2 F.1. Explicaci´on del programa Vamos ahora a presentar un m´etodo alternativo al programa anterior que consistir´a en construir la matriz del Hamiltoniano directamente en el espacio de los autovectores con un n´umero “n” de excitaciones. Primero definimos como antes los par´ametros: dim1, n, sitios,ω1, ω2, Ω, J1 y J2. Ω = −2.; (*Ω ≡∆Ω*) ω1=6.; ω2=2.∗Ω + ω1; J1 = −5; J2 = p2∗J1 ∗(J1 −Ω); Tambi´en: sitios = 3; dim = 2; dim2 = dim ∗1.0; n= dim2 −1; Definimos la base del subespacio con “n” excitaciones: Basis = SparseArray Table 1,0∗IntegerDigits[m, dim,sitios],m, 0,(dim)sitios −1; Res = SparseArray Table If[Sum[Basis[[j]][[i]],{i, 1,sitios}] == n, j −1,0,0,0],j, 1,(dim)sitios; W= Select[Res,#6= 0,0&]; Vec = Table If[Sum[Basis[[j]][[i]],{i, 1,sitios}] == n, Basis[[j]][[]],0,0],j, 1,(dim)sitios; Vec2 = Select[Vec,#6= Cero&]; L= Length[Vec2]; Definimos de forma an´aloga al caso anterior en Basis la base del espacio completo, en Wdefinimos una lista de n´umeros que corresponden a todos los autoestados de un determinado n´umero de excitaciones por ejemplo supongamos que en Basis tenemos la base completa de una excitaci´on como m´aximo por sitios y un conjunto de 3 sitios tendr´ıamos en este caso:             000 001 010 011 100 101 110 111             1En donde dim2, juega el mismo papel que dim, solo que lo definimos como un real y no como un entero para agilizar los c´alculos. 59 Como se ve las filas son los estados: |0,0,0i,|0,0,1i,|0,1,0i, ..., notemos que estos los podemos interpretar como n´umeros en base 2: 000,001,010, ..., (aqu´ı hay como m´aximo una excitaci´on por sitios pero si hubiera “n” excitaciones ser´ıa una lista de n´umeros en base “n” desde el 000 hasta el nnn) que para agilizar los c´alculos podemos usar directamente sus n´umeros decimales: {0,1,2,3,4,5,6,7}. Con ello en Wes la lista de n´umeros decimales2 de esa lista que corresponden a un determinado estado con “n” excitaciones en nuestro caso 1 excitaci´on. En nuestro caso es:   1 2 4 ↔  001 010 100  En cuanto a la variable Vec2 corresponde con el conjunto de estados con “n” excitaciones y no sus correspondientes n´umeros decimales que es W. En nuestro caso es la parte derecha de la anterior f´ormula es decir: {|0,0,1i,|0,1,0i,|1,0,0i}. A continuaci´on definimos los operadores creaci´on Ay destrucci´on a: A= Table (dim2)sitios−l,{l, 1,sitios}; a= Table −(dim2)sitios−l,{l, 1,sitios}; Primero recordemos que: c|ni=√n|n−1i(F.1) c†|ni=√n+ 1 |n+ 1i(F.2) Es decir que por ejemplo si: c† 1|1,1,0i=√2|2,1,0i´o c2|1,1,0i=√1|1,0,0i. Con Ayasolo modificaremos el estado | i el coeficiente √lo a˜nadiremos a parte, podemos as´ı darnos cuenta que el operador c† 1en nuestro caso equivale a sumar el n´umero en base 3: 100 al n´umero en base 3 (estado de partida) 110 obteniendo 100 + 110 = 210 as´ı mismos el caso de c2equivale a 110 −010 = 100. Ambos casos para simplificar los c´alculos y no tener que sumar o restar vectores podemos hacerlo en su equivalente decimal as´ı: 9 + 12 = 21 y 12 −3 = 9. Es decir los operadores equivalen a los n´umeros: c† 1≡9 y c2≡ −3. Por ello que las variables Ayaequivalen a las listas de n´umeros que representan todos los operadores creaci´on y destrucci´on respectivamente, as´ı en nuestro ejemplo de 3 sitios y una excitaci´on (base 2) ser´ıan: A=  4 2 1 a= −4 −2 −1 (F.3) Para continuar con el programa debemos de recordar la forma del Hamiltoniano: H= N X j=1 ωac† 2jc2j+ωbc† 2j−1c2j−1+J N X j=1 c† 2j−1c2j+1 +c† 2j+1c2j−1+J0 N X j=1 c† 2jc2j−1+c† 2jc2j+1 +c† 2j−1c2j+c† 2j+1c2j Para obtener la matriz del Hamiltoniano en la base del subespacio de “n” excitaciones lo que haremos ser´a dividirlo en la suma de operadores m´as sencillos, consideraremos cada pareja de operadores de creaci´on y destrucci´on como un ´unico operador: c† 2jc2j, c† 2j−1c2j−1, c† 2j−1c2j+1, c† 2j+1c2j−1, c† 2jc2j−1, c† 2jc2j+1, c† 2j−1c2jyc† 2j+1c2j. Para obtener la matriz del Hamiltoniano obtendremos todas las im´agenes de los elementos de la base, para ello haremos actuar los operadores anteriores sobre cada elemento de la base sumando a cada componente del vector Wlas correspondientes componentes de los vectores Aya. Como estos operadores conservan el n´umero de excitaciones la suma de estos tres n´umeros (A,ayW) ser´a un n´umero que este en la lista dada por Wy por tanto un vector de la base. Ser´a en el vector Imag-j (nombrando a los operadores anteriores j= 1, ..., 6)3en donde guardaremos las im´agenes de cada uno de los elementos de la base al actuar el operador j. Sin embargo recordemos que c3|0,1,0i= 0 y con nuestro m´etodo: 010 −001 = 001 ↔2−1 = 1. Para solucionar esto haremos uso de la variable Imagconst-j que calcular´a los factores √×√del la im´agenes del operador jhaciendo uso de la matriz Vec2 que guarda los vectores de los elementos de la base en sus filas, de manera que si se da el caso que si el m´etodo anterior arroja un n´umero que no corresponde realmente con la imagen del vector como cuando el operador destrucci´on act´ua sobre un sitio 2La variable Res es una variable intermedia entre Basis yW 3Notar que hay un conjunto de 8 operadores pero el caso de los dos primeros es especial y son matrices diagonales que podremos construir directamente, nuestro problema importante ser´a determinar las 6 ´ultimas matrices, para cada valor de j. 60 donde no hay excitaciones entonces el factor es cero y cuando construyamos la matriz y multipliquemos los elementos de la misma por estos factores obtendremos un cero, asegur´andonos as´ı que el m´etodo es correcto. Adem´as por asegurarnos todo elemento que sea negativo en el nuevo vector Imag-j se cambiar´a autom´aticamente por cero, as´ı tenemos: Imagconst1 = Table Table (Vec2[[m]][[2 ∗j−1]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j+ 1]])0,5,{m, 1, L},j, 1,sitios+1 2−1; Imag1 = Table Table If Imagconst1[[j]][[m]] <10−10,0,0,(A[[2 ∗j−1]] + a[[2 ∗j+ 1]] + W[[m]]), {m, 1, L}],j, 1,sitios+1 2−1; Imagconst2 = Table Table (Vec2[[m]][[2 ∗j+ 1]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j−1]])0,5,{m, 1, L},j, 1,sitios+1 2−1; Imag2 = Table Table If Imagconst2[[j]][[m]] <10−10,0,0,(a[[2 ∗j−1]] + A[[2 ∗j+ 1]] + W[[m]]), {m, 1, L}],j, 1,sitios+1 2−1; Imagconst3 = Table Table (Vec2[[m]][[2 ∗j]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j−1]])0,5,{m, 1, L},j, 1,sitios−1 2; Imag3 = Table Table If Imagconst3[[j]][[m]] <10−10,0,0,(A[[2 ∗j]] + a[[2 ∗j−1]] + W[[m]]),{m, 1, L},j, 1,sitios−1 2; Imagconst4 = Table Table (Vec2[[m]][[2 ∗j]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j+ 1]])0,5,{m, 1, L},j, 1,sitios−1 2; Imag4 = Table Table If Imagconst4[[j]][[m]] <10−10,0,0,(A[[2 ∗j]] + a[[2 ∗j+ 1]] + W[[m]]),{m, 1, L},j, 1,sitios−1 2; Imagconst5 = Table Table (Vec2[[m]][[2 ∗j−1]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j]])0,5,{m, 1, L},j, 1,sitios−1 2; Imag5 = Table Table If Imagconst5[[j]][[m]] <10−10,0,0,(A[[2 ∗j−1]] + a[[2 ∗j]] + W[[m]]),{m, 1, L},j, 1,sitios−1 2; Imagconst6 = Table Table (Vec2[[m]][[2 ∗j+ 1]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j]])0,5,{m, 1, L},j, 1,sitios−1 2; Imag6 = Table Table If Imagconst6[[j]][[m]] <10−10,0,0,(A[[2 ∗j+ 1]] + a[[2 ∗j]] + W[[m]]),{m, 1, L},j, 1,sitios−1 2; A continuaci´on lo que hacemos es comparando los distintos vectores Imag-j con el vector original Wpara obtener la matriz que representa a cada uno de los 6 operadores en la base de partida multiplicando cada elemento de la matriz por los factores de Imagconst-j. Esto lo haremos para cada valor de j y luego sumaremos todas las matrices en cada caso: MatImag1 = SparseArray Table Sum If[Imag1[[j]][[m]] == W[[l]],Imagconst1[[j]][[m]],0,0],j, 1,sitios+1 2−1, {l, 1, L},{m, 1, L}]]; MatImag2 = SparseArray Table Sum If[Imag2[[j]][[m]] == W[[l]],Imagconst2[[j]][[m]],0,0],j, 1,sitios+1 2−1, {l, 1, L},{m, 1, L}]]; MatImag3 = SparseArray Table Sum If[Imag3[[j]][[m]] == W[[l]],Imagconst3[[j]][[m]],0,0],j, 1,sitios−1 2, {l, 1, L},{m, 1, L}]]; MatImag4 = SparseArray Table Sum If[Imag4[[j]][[m]] == W[[l]],Imagconst4[[j]][[m]],0,0],j, 1,sitios−1 2, {l, 1, L},{m, 1, L}]]; MatImag5 = SparseArray Table Sum If[Imag5[[j]][[m]] == W[[l]],Imagconst5[[j]][[m]],0,0],j, 1,sitios−1 2, {l, 1, L},{m, 1, L}]]; MatImag6 = SparseArray Table Sum If[Imag6[[j]][[m]] == W[[l]],Imagconst6[[j]][[m]],0,0],j, 1,sitios−1 2, {l, 1, L},{m, 1, L}]]; Por ´ultimo solo nos quedan las matrices de los operadores c† 2jc2jy c† 2j−1c2j−1que son dos matrices diagonales en la que en la primera los elementos impares de la diagonal son nulos y viceversa en la segunda, haci´endolo para cada valor de jobtenemos: MatId1 = SparseArray Table Sum (If[l== m, 1,0,0,0]) ∗Vec2[[m]][[2 ∗j−1]],j, 1,sitios−1 2+ 1,{l, 1, L},{m, 1, L}; MatId2 = SparseArray Table Sum (If[l== m, 1,0,0,0]) ∗Vec2[[m]][[2 ∗j]],j, 1,sitios−1 2,{l, 1, L},{m, 1, L}; 61 Para finalizar solo nos queda definir el Hamiltoniano multiplicando y sumando todas las matrices y a continuaci´on obtener los autovalores y autovectores: H= SparseArray[ω1∗MatId2 + ω2∗MatId1 + J1 ∗(MatImag1 + MatImag2) + J2 ∗(MatImag3 + MatImag4 + MatImag5 + MatImag6)]; Eigenvalues[H] Eigenvectors[H] Si obtenemos la dimensi´on de la matriz para comparar con el m´etodo anterior para ver si realmente hemos simplificado el c´alculo. Para ello necesitamos calcular el n´umero de elementos de la base dado un “N” y un “n”. El problema se puede simplificar a el n´umero de combinaciones de meter “n” bolas (excitaciones) en “N” cajas (“sitios”), o lo que es lo mismo ordenar (N-1) “separadores” y n “bolas” es decir (N−1+n)!, pero como podemos intercambiar dos bolas o separadores entre s´ı sin cambiar de estado tenemos entonces que quitar las permutaciones de estos obteniendo: (N−1 + n)! n!(N−1)! =N−1 + n n Que si lo comparamos con (n+ 1)Nse puede demostrar que N−1 + n n<(n+ 1)Nde hecho para n= 5 y N= 7 tenemos que 7−1+5 5= 462 (5 + 1)7= 279936. Esto solo es el tama˜no del lado de la matriz cuadrada el total de n´umeros de cada matriz corresponde al cuadrado de esas cantidades. 62 F.2. Resumen del c´odigo El resumen de todo el c´odigo usado es: Ω = −2.; (*Ω ≡∆Ω*) ω1=6.; ω2=2.∗Ω + ω1; J1 = −5; J2 = p2∗J1 ∗(J1 −Ω); sitios = 3; dim = 2; dim2 = dim ∗1.0; n= dim2 −1; Basis = SparseArray Table 1,0∗IntegerDigits[m, dim,sitios],m, 0,(dim)sitios −1; Res = SparseArray Table If[Sum[Basis[[j]][[i]],{i, 1,sitios}] == n, j −1,0,0,0],j, 1,(dim)sitios; W= Select[Res,#6= 0,0&]; Vec = Table If[Sum[Basis[[j]][[i]],{i, 1,sitios}] == n, Basis[[j]][[]],0,0],j, 1,(dim)sitios; Vec2 = Select[Vec,#6= Cero&]; L= Length[Vec2]; A= Table (dim2)sitios−l,{l, 1,sitios}; a= Table −(dim2)sitios−l,{l, 1,sitios}; Imagconst1 = Table Table (Vec2[[m]][[2 ∗j−1]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j+ 1]])0,5,{m, 1, L},j, 1,sitios+1 2−1; Imag1 = Table Table If Imagconst1[[j]][[m]] <10−10,0,0,(A[[2 ∗j−1]] + a[[2 ∗j+ 1]] + W[[m]]), {m, 1, L}],j, 1,sitios+1 2−1; Imagconst2 = Table Table (Vec2[[m]][[2 ∗j+ 1]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j−1]])0,5,{m, 1, L},j, 1,sitios+1 2−1; Imag2 = Table Table If Imagconst2[[j]][[m]] <10−10,0,0,(a[[2 ∗j−1]] + A[[2 ∗j+ 1]] + W[[m]]), {m, 1, L}],j, 1,sitios+1 2−1; Imagconst3 = Table Table (Vec2[[m]][[2 ∗j]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j−1]])0,5,{m, 1, L},j, 1,sitios−1 2; Imag3 = Table Table If Imagconst3[[j]][[m]] <10−10,0,0,(A[[2 ∗j]] + a[[2 ∗j−1]] + W[[m]]),{m, 1, L},j, 1,sitios−1 2; Imagconst4 = Table Table (Vec2[[m]][[2 ∗j]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j+ 1]])0,5,{m, 1, L},j, 1,sitios−1 2; Imag4 = Table Table If Imagconst4[[j]][[m]] <10−10,0,0,(A[[2 ∗j]] + a[[2 ∗j+ 1]] + W[[m]]),{m, 1, L},j, 1,sitios−1 2; Imagconst5 = Table Table (Vec2[[m]][[2 ∗j−1]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j]])0,5,{m, 1, L},j, 1,sitios−1 2; Imag5 = Table Table If Imagconst5[[j]][[m]] <10−10,0,0,(A[[2 ∗j−1]] + a[[2 ∗j]] + W[[m]]),{m, 1, L},j, 1,sitios−1 2; Imagconst6 = Table Table (Vec2[[m]][[2 ∗j+ 1]] + 1,0)0,5∗(Vec2[[m]][[2 ∗j]])0,5,{m, 1, L},j, 1,sitios−1 2; Imag6 = Table Table If Imagconst6[[j]][[m]] <10−10,0,0,(A[[2 ∗j+ 1]] + a[[2 ∗j]] + W[[m]]),{m, 1, L},j, 1,sitios−1 2; MatImag1 = SparseArray Table Sum If[Imag1[[j]][[m]] == W[[l]],Imagconst1[[j]][[m]],0,0],j, 1,sitios+1 2−1, {l, 1, L},{m, 1, L}]]; MatImag2 = SparseArray Table Sum If[Imag2[[j]][[m]] == W[[l]],Imagconst2[[j]][[m]],0,0],j, 1,sitios+1 2−1, {l, 1, L},{m, 1, L}]]; MatImag3 = SparseArray Table Sum If[Imag3[[j]][[m]] == W[[l]],Imagconst3[[j]][[m]],0,0],j, 1,sitios−1 2, {l, 1, L},{m, 1, L}]]; 63 MatImag4 = SparseArray Table Sum If[Imag4[[j]][[m]] == W[[l]],Imagconst4[[j]][[m]],0,0],j, 1,sitios−1 2, {l, 1, L},{m, 1, L}]]; MatImag5 = SparseArray Table Sum If[Imag5[[j]][[m]] == W[[l]],Imagconst5[[j]][[m]],0,0],j, 1,sitios−1 2, {l, 1, L},{m, 1, L}]]; MatImag6 = SparseArray Table Sum If[Imag6[[j]][[m]] == W[[l]],Imagconst6[[j]][[m]],0,0],j, 1,sitios−1 2, {l, 1, L},{m, 1, L}]]; MatId1 = SparseArray Table Sum (If[l== m, 1,0,0,0]) ∗Vec2[[m]][[2 ∗j−1]],j, 1,sitios−1 2+ 1,{l, 1, L},{m, 1, L}; MatId2 = SparseArray Table Sum (If[l== m, 1,0,0,0]) ∗Vec2[[m]][[2 ∗j]],j, 1,sitios−1 2,{l, 1, L},{m, 1, L}; H= SparseArray[ω1∗MatId2+ω2∗MatId1+J1∗(MatImag1+MatImag2)+J2∗(MatImag3+MatImag4+MatImag5+ MatImag6)]; Eigenvalues[H] Eigenvectors[H] 64 Ap´endice G C´odigo M´etodo 2 - Modificaci´on para comprobar los auto estados localizados No tenemos m´as que a˜nadir el siguiente c´odigo al final del c´odigo del Ap´endice F. Identidad = Table[Table[KroneckerDelta[j, m],{j, 1,sitios}],{m, 1,sitios}]; AutoVector = Table Identidad[[2 ∗j]] −√2∗Identidad[[2 ∗j+ 1]] + Identidad[[2 ∗j+ 2]]∗0,5,j, 1,sitios−3 2; MatrixForm[AutoVector]Normalizado = Table Sum[AutoVector[[j]][[m]] ∗AutoVector[[j]][[m]],{m, 1,sitios}],j, 1,sitios−3 2; MatrixForm[Normalizado]AutoVectorNorma = Table Table[AutoVector[[j]][[m]]/Normalizado[[j]],{m, 1,sitios}],j, 1,sitios−3 2; MatrixForm[AutoVectorNorma]Imag = Table Table[Sum[H[[m]][[l]] ∗AutoVector[[j]][[l]],{l, 1,sitios}],{m, 1,sitios}],j, 1,sitios−3 2; MatrixForm[Imag]AutoValTeo = Table Sum[Imag[[j]][[m]] ∗AutoVector[[j]][[m]]/Normalizado[[j]],{m, 1,sitios}],j, 1,sitios−3 2; MatrixForm[AutoValTeo] 65 66 Ap´endice H C´odigo M´etodo 2 - Modificaci´on para comprobar los auto estados localizados, frecuencias distintas Deberemos a˜nadir lo siguiente al c´odigo expuesto en el Ap´endice F. Primero cambiaremos la condici´on inicial para J1: J1 = −1+√2Ω+2Ω2−√1+2√2Ω+4Ω2+4√2Ω3+4Ω4 2√2+3Ω Segundo a˜nadimos el siguiente c´odigo al final: Identidad = Table[Table[KroneckerDelta[j, m],{j, 1,sitios}],{m, 1,sitios}]; AutoVector = Table hqJ1 J1−Ω∗Identidad[[2 ∗j]] −√2∗Identidad[[2 ∗j+ 1]] +qJ1 J1−Ω∗Identidad[[2 ∗j+ 2]]∗0,5,j, 1,sitios−3 2i; MatrixForm[AutoVector]Normalizado = Table Sum[AutoVector[[j]][[m]] ∗AutoVector[[j]][[m]],{m, 1,sitios}],j, 1,sitios−3 2; MatrixForm[Normalizado]AutoVectorNorma = Table Table[AutoVector[[j]][[m]]/Normalizado[[j]],{m, 1,sitios}],j, 1,sitios−3 2; MatrixForm[AutoVectorNorma]Imag = Table Table[Sum[H[[m]][[l]] ∗AutoVector[[j]][[l]],{l, 1,sitios}],{m, 1,sitios}],j, 1,sitios−3 2; MatrixForm[Imag]AutoValTeo = Table Sum[Imag[[j]][[m]] ∗AutoVector[[j]][[m]]/Normalizado[[j]],{m, 1,sitios}],j, 1,sitios−3 2; MatrixForm[AutoValTeo] 67