Fenomenolog´ıa Eduardo Menendez July 21, 2005

2

Contents 1 Elementos de f´ısica de s´ olidos 1.1 Elementos de cristalograf´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Ejemplos. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Simetr´ıa . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Fuerzas interat´omicas y tipos de s´olidos . . . . . . . . . . . . . . . . . . . . . . 1.2.1 “Gases nobles” . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Cristales i´onicos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.3 Cristales covalentes . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.4 Metales . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Elasticidad . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.1 Simetr´ıas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.2 Deformaciones homog´eneas . . . . . . . . . . . . . . . . . . . . . . . . . 1.3.3 Ondas el´asticas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.4 Modelos discretos de las vibraciones de la red cristalina . . . . . . . . . . . . . 1.4.1 Modos normales de una cadena lineal monoat´omica . . . . . . . . . . . 1.4.2 Redes 1-D con bases . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5 Ecuaci´on de estado . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.1 Ecuaci´on de Murnaghan . . . . . . . . . . . . . . . . . . . . . . . . . . 1.5.2 Ecuaci´on de Birch-Murnaghan . . . . . . . . . . . . . . . . . . . . . . . 1.6 Densidad de Estados . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.6.1 Particula libre en una caja unidimensional . . . . . . . . . . . . . . . . 1.6.2 Part´ıcula libre con condiciones de frontera periodicas en una dimensi´on 1.6.3 Particula libre en dos dimensiones . . . . . . . . . . . . . . . . . . . . . 1.6.4 Particula libre en tres dimensiones . . . . . . . . . . . . . . . . . . . . 1.6.5 Cadena lineal monoat´omica . . . . . . . . . . . . . . . . . . . . . . . . 1.6.6 Vibraciones de las redes tridimensionales . . . . . . . . . . . . . . . . . 1.6.7 Modelos del calor especifico . . . . . . . . . . . . . . . . . . . . . . . . 1.7 Despu´es de la aproximaci´on arm´onica . . . . . . . . . . . . . . . . . . . . . . . 1.7.1 Efectos cu´anticos . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7.2 Pasos de una simulaci´on . . . . . . . . . . . . . . . . . . . . . . . . . . 1.7.3 C´alculo de promedios y errores estad´ısticos . . . . . . . . . . . . . . . .

5 5 9 10 15 16 19 19 20 21 23 26 28 30 30 32 35 35 36 37 37 38 39 39 40 40 46 48 52 53 55

2 Elementos de ´ optica no lineal 61 2.1 Efecto electro´optico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 2.2 Mezcla de dos ondas . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 3

4

CONTENTS

2.3

2.2.1 Phase matching en problemas unidimensionales 2.2.2 Aplicaciones . . . . . . . . . . . . . . . . . . . . Fen´omenos de tercer orden . . . . . . . . . . . . . . . . 2.3.1 Efecto electro´optico . . . . . . . . . . . . . . . . 2.3.2 Tercer arm´onico y efecto Kerr ´optico . . . . . . 2.3.3 Autoenfoque . . . . . . . . . . . . . . . . . . . .

3 Elementos de biof´ısica 3.1 Biomec´anica . . . . . . . . . . . . . . . . . . . . . . . 3.1.1 Fuerzas que act´ uan sobre el f´emur y la cadera 3.1.2 Fuerzas en equlibrio sobre un pie . . . . . . . 3.1.3 Efecto de un bast´on . . . . . . . . . . . . . . 3.2 Fuerzas que act´ uan sobre las v´ertebras lumbares . . .

. . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . . . . . . .

. . . . . .

66 66 67 67 67 68

. . . . .

71 71 71 72 74 76

4 F´ısica de part´ıculas elementales

81

5 F´ısica del n´ ucleo at´ omico

83

6 Mapas discretos 85 6.1 El mapa log´ıstico . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85 7 Caos en ecuaciones diferenciales ordinarias 91 7.0.1 Corte de Poincar´e . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93 7.0.2 Mapa estrobosc´opico . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94

Chapter 1 Elementos de f´ısica de s´ olidos 1.1

Elementos de cristalograf´ıa

Los materiales s´olidos se clasifican de diversas formas, seg´ un determinados par´ametros, por ejemplo: conductores y aislantes, tranasparentes y opacos, ionicos o metalicos o covalentes, cristalinos o amorfos. Comenzaremos con esta clasificaci´on, que responde a la forma en que se ordenan los ´atomos. Un cristal consiste en atomos dispuestos segun un patron que se repite periodicamente. La mayoria de los materiales son cristales. Se pueden subdividir en monocristalinos y policristalino. En un monocristal el patr´on que se repite periodicamente es unico y se extiende a lo largo de dimensiones macroscopicas. Ejemplo son los diamantes y otras piedras preciosas, los cristales de sal, los osciladores de un reloj de cuarzo. Un material policristalino se compone de una multitud de monocristales microscopicos (llamados granos) que estan orientados en todas las direcciones posibles. Tanto monocristales como policristales tienen la caracteristica de presentar un orden casi perfecto en volumenes grandes comparados con las dimensiones atomicas. Un material amorfo es aquel que no es cristalino. Ejemplo de los materiales amorfos son los vidrios. En estos no existe un patron que se repita periodicamente. La periodicidad tiene consecuencias importantes. En los materiales cristalinos ocurren una serie de fenomenos que no se dan en los amorfos, por ejemplo, transiciones de fases bien definidas, conductividad electrica, bandas de energia y difraccion de rayos X con maximos en angulos bien determinados. La difraccion de rayos X constituye la tecnica experimental mas poderosa para determinar si un material es cristalino o amorfo. Veamos como describir un cristal de forma precisa. Red de Bravais. Rs una coleccion de puntos que llena el espacio, tales que cada punto tiene exactamente el mismo entorno. Vectores de traslacion. Es cualquier vector que conecta dos puntos de una red de Bravais. Por definicion, una red de Bravais es infinita. Si no lo fuese, los puntos de la superficie no tendrian el mismo entorno que los del interior. Si a una red de Bravais se traslada segun un vector de traslacion, la red queda invariante. Esto tiene consecuencias profundas, que se estudian mediante la teoria de grupos. 5

6

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

Figure 1.1: Dos posibilidades de escoger los vectores primitivos en una red de Bravais bidimensional. Los dos conjuntos de vectores primitivos se obtienen uno del otro. ˜t1 = 2t1 − t2 y ˜t2 = t2 − t1 . Vectores primitivos Vectores de traslacion ~a1 , ~a2 , ~a3 , tales que todos los puntos de la ~ = n1~a1 + n2~a2 + n3~a3 , donde red de Bravais se pueden conectar por vectores de traslaci´on R n1 , n2 , n3 son numeros enteros. Celda unitaria es el paralelepipedo generado por tres vectores de la red no coplanares. Si estos tres vectores son primitivos, entonces se denomina celda primitiva. En una red bidimensional, las celdas unitarias son paralelogramos generado por dos vectores de la red no colineales. Las celdas primitiva cumple que, trasladada segun todos los vectores de la red de Bravais, llena todo el espacio sin traslaparse ni dejar vacio. Los vectores primitivos no son unicos. Considere la siguiente red bidimensional de la Figura 1.1. Hay infinitas formas de escoger los vectores pimitivos, aunque unas son mas comodas que otras. La caracteristica que los une es que todas las celda primitivas generada por ellos tiene la misma area. En una red tridimensional, todas las celdas primitivas tienen el mismo volumen. El volumen de una celda definita por vectores ~a1 , ~a2 , ~a3 es igual a Ω = ~a1 · ~a2 × ~a3 .

(1.1)

Toda celda primitiva es una celda unitaria. Una celda unitaria de volumen minimo, es una celda primitiva. Notemos que a una celda primitiva le corresponde exactamente un punto de la red. En la Fig. 1.1 podemos notar que cada celda tiene puntos en los vertices. Cada uno de estos puntos esta compartido por cuatro celdas (la celda central y tres vecinas), de modo que cada uno contribuye un cuarto a la celda central. Sumando se obtiene un punto. Una celda unitaria no primitiva contiene m´as de un punto de la red. En ocasiones la celda primitiva no tiene toda la simetria de la red. En estos casos se acostumbra utilizar celdas que tiene todas las propiedades de simetria de la red. Generalmente esta celda tiene mayor volumen y se le llama celda convencional.

1.1. ELEMENTOS DE CRISTALOGRAF´IA

7

Seg´ un hemos visto, especificando tres vectores primitivos se define completamente una red de Bravais, y a su vez se define la celda primitiva. Es usual en cristalografia describir la celda primitiva mediante 6 par´amtros: a, b, c, α, β, γ, donde a, b y c son los lados de la celda primitiva, α es el ´angulo formado por los lados b y c, β por los lados a y c, y γ por los lados a y b. Esto se ilustra en la Fig. 1.2. Notese que definir tres vectores requiere 9

Figure 1.2: numeros, mientras que solo se necesitan 6 para definir la celda unitaria. Los 3 parametros que sobran en la definicion de los vectores definen la orientacion de la celda respecto al sistema de coordenadas. En la practica, se escoge el sistema de coordenadas mas comodo posible. Muchos programas de cristalografia o de calculos atomisticos definen el vector ~a1 = (a, 0, 0) paralelo al eje OX (1 parametro), ~a2 = (b cos γ, b sin γ, 0) (2 parametros) y el vector ~a3 en una direcci´on general (3 parametros). Ejercicio. Encuentre una formula para las componentes de este vector!.

Figure 1.3: Red hexagonal que no corresponde a una red de Bravais, pues los puntos blancos no tienen el mismo entorno que los puntos negros. Tambi´en se muestra la red de Bravais asociada. Notemos que no toda red es de Bravais, lo cual se ilustra el la Fig. 1.3. Sin embargo, se le puede asociar una “red asociada” a los puntos negros o a los puntos blancos. La figura anterior se corresponde al ordenamiento de los atomos de carbono en ciertos planos

8

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

del grafito. La estructura tridimensional del grafito se obtiene superponiendo estos planos en direccion perpendicular a la figura. La estructura bidimensional ideal ilustrada en la figura se denomina grafeno, t´ermino muy utilizado en la literatura de nanotubos de carbono. Este caso demuestra que no hay cristales (la mayoria) que no se pueden describir solamene con una red de Bravais. Esto hace necesario el concepto de base. Base Conjunto de atomos que se asocia identicamente a cada punto de la red. Generalmente la base se especifica mediante las coordenadas de los ´atomos respecto a un punto de la red, que se toma como origen de coordenadas, y respecto a ciertos ejes de coordenadas. Si los vectores primitivos se especifican mediante sus componentes, entonces las posiciones at´omicas se especifican mediante sus coordenadas cartesianas, respecto al mismo sistema de ejes usado para definir los vectores. Por otra parte, si la celda (primitiva o convencional) se define mediante los par´ametros a, b, c, α, β, γ, entonces las posiciones at´omicas se especifican en coordenadas fraccionarias x, y, z, de modo que la posici´on de un at´omo es ~r = x~a1 + y~a2 + z~a3 ,

(0 ≤ x, y, z < 1).

(1.2)

La red del grafeno es un claro ejemplo de red con base, siendo la red, por ejemplo, los puntos negros en la Fig. 1.3, y la base el conjunto de dos atomos ubicados en un punto negro y en el punto blanco superior. Para aclarar el concepto de base, veamos la simpatica Figura 1.4. Las cruces × son puntos de una red, mientras que la flor constituye la base asociada a cada punto de la red. N´otese que los puntos de la red son puntos matem´aticos, no es necesario que esten ocupados por ´atomos o por ninguna parte de la base.

Figure 1.4: Ejemplo de una red en la cual el los simbolos × indican los puntos de la red, y la flor es la base.

1.1. ELEMENTOS DE CRISTALOGRAF´IA

1.1.1

9

Ejemplos.

Red cubica simple Los vectores primitivos ~a1 = aˆ x ~a2 = aˆ y ~a3 = aˆ z

(1.3)

donde a es el par´ametro de red. La celda primitiva es un cubo. El Polonio presenta esta estructura. Red c´ ubica de cara centrada (fcc).

Figure 1.5: Red c´ ubica de cara centrada (fcc) con par´ametro de red a. Derecha: Celda convencional y celda primitiva. Los vectores primitivos son a y + zˆ) ~a1 = (ˆ 2 a ~a2 = (ˆ z + xˆ) 2 a ~a3 = (ˆ x + yˆ) 2

(1.4)

donde a es el par´ametro de red. El Cobre y el Arg´on presentan esta estructura con par´ametros de red aCu = 3.61 [˚ A] y aAr = 5.26 [˚ A] a (4.2 [K]) respectivamente. La Fig. 1.5 muestra la celda convencional de la red fcc. Contemos el numero de puntos de la red que contiene esta celda. Hay 6 puntos en el centro de las caras, que se comparten entre dos celdas adyacentes, por tanto, contamos 3 puntos. En los vertices del cubo hay 8 puntos, cada uno de los cuales se comparten en 8 celdas, por tanto, suman 1 punto. Asi en total corresponden 4 puntos de la red. Esto indica que el volumen de la celda convencional a3 es 4 veces el volumen de la celda primitiva, definida por los vectores (1.4). La Figura 1.5 muestra la relacion entre la celda convencional y la celda primitiva. Notese que la celda convencional tiene la simetria cubica de la red, pero no asi´ı la celda primitiva. Ejercicio Demuestre usando la formula (1.1) que el volumen de la celda primitiva es a3 /4. Definimos el n´ umero de coordinaci´ on como el n´ umero de vecinos m´as cercanos. Para los ejemplo tenemos:

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

10

Figure 1.6: Ilustracion de la simetria de un cubo. • C´ ubica simple seis vecinos. • C´ ubica de cuerpo centrado (bcc) ocho vecinos. • C´ ubica centrada en las caras (fcc) doce vecinos.

1.1.2

Simetr´ıa

Una operacion de simetr´ıa es aquella que que aplicada a un sistema lo deja invariante. ~ = n1~a1 + n2~a2 + n3~a3 (n1 , n2 , n3 enteros) Como ejemplo, ya hemos visto que las traslaciones R dejan invariantes a una red de Bravais. El estudio de la simetria de los cristales facilita la clasificacion de estos y el calculo de muchas magnitudes. Toda red de Bravais presenta simetria de traslacion. Adicionalmente puede presentar otras operaciones. Consideremos la red fcc, que hemos visto en los ejemplos. La celda convencional tiene todas las propiedades de simetria de un cubo. V´ease la Figura 1.6 1. Reflexion respecto a un planos paralelos a las caras (σ). 2. Reflexion respecto a planos que pasan por las diagonales de las caras. 3. Rotaciones en multiplos de 90◦ alrededor de los ejes C4 . 4. Rotaciones en multiplos de 120◦ alrededor de los ejes S6 , que pasan por las diagonales. 5. Roto-reflexiones en el eje S6 . Esto es una rotacion de 60◦ , seguida de una reflexion en el plano perpendicular al eje S6 y que pasa por el centro del cubo. 6. Rotaciones en multiplos de 180◦ alrededor de los ejes C2 . 7. Inversion respecto al centro del cubo. 8. Combinaciones de todas las operaciones anteriores. No pretendemos hacer un estudio detallado de las propiedades de simetria, sino informar de su existencias y de los conceptos basicos. El conjunto de propiedades de simetria de un sistema forma una estructura matematica llamada grupo, que cumple unos pocos axiomas

1.1. ELEMENTOS DE CRISTALOGRAF´IA

11

1. La combinacion de dos operaciones de simetria es una operacion de simetria 2. Existe la operacion identidad (no hacer nada). 3. Para toda operacion de simetria existe la operacion inversa (ejemplo, rotaciones de 90◦ y -90◦ .) A partir de las propiedades anteriores y estableciendo una relacion con el algebra de matrices se ha desarrollado la llamada teoria de grupos, que es la teoria matematica de la simetria. De los postulados anteriores puede razonarse que el conjunto de traslaciones de una red de Bravais forma un grupo, el Grupo de Traslaciones de la Red de Bravais. El conjunto de transformaciones de simetria de un sistema (por ejemplo, un cubo, una red), que incluye rotaciones, reflexiones y roto-reflexiones respecto a ejes que pasan por un punto comun, ademas de la inversion respeto al mismo punto, se denomina grupo puntual. La red fcc tiene el mismo grupo puntual de simetria que el cubo, que en este caso particular se llama Oh . La celda convencional de una red, por su definicion, tiene el mismo grupo puntual que la red de Bravais correspondiente. Si se combina una traslacion con una operacion del grupo puntual, tambien se deja invariante la red de Bravais. Al conjunto de operaciones de simetria que incluyen traslaciones y operaciones puntuales, se le llama grupo espacial de la red de Bravais. Para las redes de Bravais tridimensionales existen solamente siete grupos puntuales posibles y 14 grupos espaciales. Obviamente, varios grupos espaciales comportan el mismo grupo puntual. Esto permite clasificar todos los cristales en siete sistemas cristalinos (seg´ un el grupo puntual) y en 14 redes de Bravais (segun el grupo espacial). La Figura 1.7 ilustra los siete sistemas cristalinos 1. Cubico. a = b = c, α = β = γ = 90◦ . Tiene tres posibles redes, que se denotan con letras mayusculas seg´ un el tipo de celda convencional y los puntos de la red que esta contiene: la cubica simple (P, de Primitive), cubica centraba en el cuerpo (I, de Inner), cubica centrada en las caras (F, de Faces). Tambien se conocen por sus siglas en ingles: sc (simple cubic), bcc (body centered cubic) y fcc (face centered cubic). 2. Tetragonal a = b 6= c, α = β = γ = 90◦ . Tiene dos posibles redes, simple (P) y centrada (I). 3. Ortorrombico a 6= b 6= c, α = β = γ = 90◦ . Tiene cuatro posibles redes: simple (P), centrada en el cuerpo (I), centrada en las bases (C, de Center), centrada en las caras (F). 4. Monoclinico. a 6= b 6= c, β o γ 6= 90◦ . Tiene dos posibles redes: simple (P) y centrada en las bases (C). 5. Triclinico. Solo tiene redes simples y celdas primitivas (P). 6. Hexagonal. a = b, α = β = 90◦ , γ = 120◦ . Solo hay un tipo de red (P). La celda convencional no es un hexagono, aunque se ilustre asi en la figura. Recuerde que siempre la celda convencional es un paralelepipedo.

12

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS 7. Trigonal. a = b = c, α = β = γ < 120◦ . Para este sistema siempre es posible escoger alternativamente una celda del tipo hexagonal. Si esta es primitiva, entonces es equivalente al sistema hexagonal. Si la celda hexagonal no es primitiva, entonces la celda primitiva es romboedrica y se designa R.

Figure 1.7: Las 14 redes de Bravais, representadas por sus celdas convencionales. En 1842, M. L. Frankenheim, determino erroneamente que existen 15 tipos de redes. A.

1.1. ELEMENTOS DE CRISTALOGRAF´IA

13

Figure 1.8: Ilustracion de la simetria de la molecula de agua. Bravais corrigio el error en 1845 y por eso las redes llevan su apellido.1 La Figura 1.8 representa la molecula de agua. Las propiedades de simetria que esta tiene son: identidad, rotacion de 180◦ alrededor del eje central, reflexion respecto al plano de la molecula y respecto al plano perpendicular, y todas las combinaciones de estas. El conjunto de estas propiedades forma el grupo de simetria denominado C2v . ¿Qu´e ocurre si en cada punto de la red de Bravais fcc se pone una molecula de agua? ¿Que simetria tiene ese hipotetico cristal? Evidentemente, las operaciones de simetria son las combinaciones de las traslaciones de la red m´as aquellas operaciones puntuales que son comunes a la red fcc y a la molecula de agua. en algunos casos. La simetria de este cristal es mas baja que la simetria de la red de Bravais, y el conjunto de operaciones de simetria que dejan invariante el cristal, se llama grupo espacial del cristal o grupo espacial cristalografico. Esto implica que los subgrupos2 de los 14 grupos espaciales de las redes de Bravais, tambien son posibles para las estructuras cristalinas. Estos son todos los grupos de simetria que puede tener una red de Bravais con base y se denominan grupos espaciales cristalograficos. Existe un total de 230 grupos espaciales cristalograficos (comparese con los 14 grupos espaciales de las redes sin base). El subconjunto de operaciones de simetria puntuales (o sea, las operaciones que no involucran traslacion) de un grupo espacial cristalografico se llama grupo puntual cristalografico, de los cuales existen 32 (comparese con los 7 grupos puntuales de las redes de Bravais). No nos detendremos a explicar la notacion utilizada para designar los grupos espaciales. Basta conocer que existen 230 y que estos especifican totalmente la simetria de un cristal. En las bases de datos de cristalografia, y en muchos programas de simulacion, la informacion del grupo espacial es esencial. Existen tres formas de designarlos 1. Numerica. Simplemente un numero entre 1 y 230, segun un orden establecido en un libro llamado Tablas Internacionales de Cristalografia. 2. Notacion de Schoenflies. Ejemplo, el grupo 14 es C52h 1

Con justicia, deberian llamarse redes de Frankenheim-Bravais. En el siglo XX, A. A. Abrikosov (premio Nobel 2003) predijo las redes de vortices en los superconductores de alta temperatura. Por errores numericos de calculo, predijo que las redes deberian ser cuadradas, cuando en realidad son hexagonales. Sin embargo se conocen como redes de Abrikosov (y desconozco quien enmendo el error). 2 Un subgrupo es un subconjunto de un grupo que tambien es un grupo, o sea, que cumple los tres axiomas de grupo.

14

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS 3. Notacion de Herman-Maughin. Ejemplo, el grupo 14 es P 1 1 21 /b 1, o abreviadamente P21 /b.

El trabajo de deducir las propiedades de todos los grupos es monumental y se encuentra compilado en las Tablas Internacionales de Cristalografia, las que constituyen una herramienta fundamental para el trabajo de los cristalografos. Un fragmento de esta tabla se muestra en la Figura 1.9

Figure 1.9: Fragmento de las Tablas Internacionales de Cristalografia Una nota respecto a la informacion contenida en estas tablas. En la parte izquierda puede verse una columna con el encabezamiento “Number of positions, Wyckoff notation, and point

´ ´ 1.2. FUERZAS INTERATOMICAS Y TIPOS DE SOLIDOS

15

symmetry”. Estos datos indican las posiciones en que pueden encontrarse atomos dentro de la celda convencional (en este caso es primitiva, indicado por la letra P del nombre del grupo). A su derecha se indican las coordenadas de estas posiciones. Todas las coordenadas asociadas a un tipo de posicion son equivalentes por simetria. En una estructura cristalina, una posicion puede estar ocupada por un atomo o vacia, pero si esta ocupada, estan ocupadas todas las posiciones equivalentes. Los programas que se usan en cristalgrafia y en simulaciones, suelen aceptar como entrada solo una de las posiciones equivalentes, siendo generadas las demas por el programa. En las bases de datos de cristalografia, y en la revista Acta Crystallographica, solo se reporta una de las posiciones equivalentes, siendo responsabilidad del usuario generar las demas.

1.2

Fuerzas interat´ omicas y tipos de s´ olidos

El problema del c´alculo de los estados estacionarios en una porci´on de un material s´olido involucra, en principio, la resoluci´on de la ecuaci´on de Schr¨odinger del cristal, en la que deben tenerse en cuenta todas las interacciones presentes y las coordenadas de todos los electrones {ri } y todos los n´ ucleos {Rα }. HΨ({Rα }, {ri }) = EΨ({Rα }, {ri }).

(1.5)

El Hamiltoniano del sistema contiene la energ´ıa cin´etica de los iones, la de los electrones de valencia y las energ´ıas de interacci´on ion-ion, ion-electr´on y electr´on-electr´on. H = Telec + Tion + Vion−elec + Vion−ion + Velec−elec X ~2 X Zα e2 ~2 2 X 2 = − ∇ + − ∇ − 2m i 2Mα α 4πε0 |ri − Rα | α i i,α +

1 X Zα Zβ e2 1X e2 + . 2 α,β6=α 4πε0 |Rβ − Rα | 2 i,j6=i 4πε0 |ri − rj |

(1.6)

donde los terminos son, por orden, los operadores de energ´ıa cin´etica de los electrones, energ´ıa cin´etica de los n´ ucleos, energ´ıa de interacci´on electr´on-n´ ucleo, energ´ıa de repulsi´on de los n´ ucleos, y energ´ıa de interacci´on entre los electones. El enorme n´ umero de part´ıculas que contienen los cristales hace imposible la resoluci´on directa de dicha ecuaci´on. Por esto se hace necesario aplicar la estrategia de divide et impera, por medio de aproximaciones. Una primera aproximaci´on es sugerida por el hecho de que la masa del electron es 1800 veces menor que la del proton, por lo cual la dinamica del movimiento se pude separar en variables rapidas (electrones) y variables lentas (nucleos). En esto consiste la llamada aproximaci´ on adiab´ atica o aproximaci´ on de Born-Oppenheimer. Como los electrones se mueven mucho m´as r´apido, puede considerarse que en cada posicion de los n´ ucleos {Rα }, los electrones se encuentran en un estado estacionario (casi siempre el estado basico) correspondiente al operador hamiltoniano Helec = Telec + Vion−elec + Velec−elec ,

(1.7)

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

16

donde las posiciones de los n´ ucles aparecen solamente como parametros y las variables de la ecuacion de Schrodinger son las coordenadas electronicas. Helec Ψn ({ri }, {Rα }) = En ({Rα })Ψn ({ri }, {Rα }).

(1.8)

Las energ´ıas En ({Rα }) actuan como un campo de fuerzas para los nucleos, obedecen al Hamiltoniano Hion = Tion + Vion−ion + En ({Rα }). (1.9) As´ı el efecto de los electrones sobre el conjunto de iones se manifiesta a trav´es de una energ´ıa potencial efectiva. As´ı, el movimiento de los electrones se desacopla del de los iones. Si bien el caracter cuantico de los electrones es fundamental para describir su comportamiento, los nucleos, por ser mas pesados y tener una longitud de onda mucho menor (h/p), pueden considerarse clasicos en la mayoria de los casos practicos (efectos cuanticos finos son observables en el hidrogeno y en los demas elementos a temperaturas cercanas al cero absoluto, notablemente el He superfluido). Cl´asicamente, los n´ ucleos obedecen las ecuaciones de Newton ¨ α = −∇α [En ({Rα }) + Vion−ion ({Rα })] = Fα . Mα R (1.10) El c´alculo de las energ´ıas electr´onicas En ({Rα }) y sus derivadas las fuerzas Fα puede ser computado a partir de la ecuacion cu´antica (1.8) o su equivalente en la llamada Teoria del Funcional de la Densidad (DFT, siglas en ingl´es), las cuales pueden ser resueltas num´ericamente en sistemas de pocos atomos3 . Dado el largo tiempo de c´alculo y los pocos ´atomos que se pueden tratar, es conveniente usar expresiones anal´ıticas que aproximen las verdaderas fuerzas. Con el trabajo de medio siglo se han determinado expresiones adaptadas para metales, s´olidos covalentes, gases nobles, l´ıquidos, etc.

1.2.1

“Gases nobles”

A este grupo corresponden los solidos formados al enfriar los llamados gases nobles (Ne, Ar, Kr, Xe). En estado atomico estos elementos tienen capas electronicas totalmente llenas y por eso tienen muy poca posibilidad de formar enlaces. Los atomos de estos elementos ineractuan mediante fuerzas de van der Waals, tambien llamadas fuerzas de dipolos fluctuantes. Considerense dos atomos (1 y 2) separados una distancia r. El dipolo instantaneo en el atomo 1 p~1 , crea un campo electrico proporcional a E ∼ p1 /r3 , que induce dipole en 2 proporcional al campo, p2 = αE ∼ αp1 /r3 . La energia de interaccion electrostatica entre los dos dipolods es del orden de αp21 p2 p1 (1.11) Vdip−dip ∼ − 3 ∼ − 6 . r r Aunque promediado en el tiempo h~p1 i = 0, la energia de interaccion es proporcional a hp21 i 6= 0. La mecanica cuantica es necesaria para calcular las constantes de proporcionalidad de forma precisa, pero el razonamiento anterior predice correctamente las dependencia proporcional a 1/r6 de la energia de interaccion. 3 El record hace unos pocos a˜ nnos era de 64000 atomos, con muchas aproximaciones y un poderoso computador. Rutinariamente se calculan sistemas de menos de 100 ´atomos con la DFT y menos de 10 en metodos m´as exactos.

´ ´ 1.2. FUERZAS INTERATOMICAS Y TIPOS DE SOLIDOS Table 1.1: Ne ² (meV) 3.1 σ (˚ A) 2.74

17

Valores de los parametros de Lennard-Jones para los gases nobles. Ar Kr Xe 10.4 14.0 20.0 3.40 3.65 3.98

Cuando los atomos estan muy proximos, de modo que las nubes electronicas entran en contacto y se interpenetran, aparecen fuerzas de repulsi´on que sobrepasan a las fuerzas de van der Waals e impiden que los atomos se junten. Las fuerzas de repulsion se explican solo desde la mecanica cuantica, calulando En ({Rα }) en la Ec. (1.8). Hay varias formas de ajustar este potencial de repulsi´on, por ejemplo, B/r12 . Con esto, el potencial de interacci´on queda escrito en la forma ·³ ´ ¸ B A σ 12 ³ σ ´6 σ = (B/A)1/6 . (1.12) − , φ(r) = 12 − 6 = 4² ² = A2 /4B r r r r La expresi´on (1.12) es conocida como potencial de Lennard-Jonnes. La potencia 12 en el t´ermino repulsivo es no se deduce, sino que se usa por simplicidad analitica. Este potencial permite reproducir las propiedades estructurales y termodinamicas de los gases nobles tanto en estado gaseoso como solido. Consideremos algunas de estas propiedades. ~ los sitios Los solidos de gases nobles forman una red fcc con bases monoatomicas. Sean R de la red, entonces la energ´ıa de interaccion de un atomo con el resto de la red es X φ(R). (1.13) V1 = ~ =0 R6

Si multiplicamos por N , el numero de atomos del cristal4 , obtenemos el doble de la energ´ıa potencial del cristal, pues se ha contado doblemente la interacci´on. as´ı, la energ´ıa por atomo se obtiene dividiendo por N/2 1X u= φ(R). (1.14) 2 ~ =0 R6

~ como un numero adimensional αR ~ Es conveniente escribir la longitud del vector de la red R multiplicado por la distancia a los primeros vecinos r. Con esto, la ecuacion (1.14) se expresa en la forma · ³ σ ´6 ¸ ³ σ ´12 − A6 , (1.15) u = 2² A12 ( r r donde An =

X ~ =0 R6

1 . ~ n α(R)

(1.16)

4 en un red infinita, N es infinito, pero para obtener la energia de enlace por atomos, dividimos por N y se cancelan los infinitos. En realidad aqui los que se esta haciendo es despreciar los efectos de superficie. Un cristal real siempre tendra N finito. El numero de atomos en la superficie es proporcional a N 2/3 y para estos la energia potencial es distinta. Como lim N 2/3 /N → 0, cuando N representa una cantidad macroscopica de atomos la energia por atomo tiende a la energia de los atomos interiores.

18

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

Table 1.2: Distancia a primeros vecinos y energ´ıa de cohesi´on de los solidos de gases nobles. Ref. Ashcroft-Mermin, Cap. 20, pag 401. Ne Ar Kr Xe ˚ r0 (A) Experimento 3.13 3.75 3.99 4.33 1.09σ 2.99 3.71 3.98 4.34 u0 (eV/atomo) Experimento -0.02 -0.08 -0.11 -0.17 −8.6² -0.027 -0.089 -0.120 -0.172 10 2 B (10 dina/cm ) Experimento 1.1 2.7 3.5 3.6 3 75²/σ 1.81 3.18 3.46 3.81 Un aproximaci´on grosera se obtiene despreciando las energ´ıas de interacci´on entre los atomos mas alejados que la distancia de primeros vecinos r. en el caso de la red fcc, el numero ~ An = 12. Si usamos los de primeros vecinos es 12 y por tanto, dada la definici´on de αR, vectores primitivos dados por la ecuacion (1.4) y hacemos la suma num´erica sobre toda la red, se obtiene A6 = 14.45 y A12 = 12.13. Ejercicio Eval´ ue las constantes A6 y A12 para las redes fcc, bcc, y sc. La constante de la red y la densidad de equilibrio se obtienen f´acilmente minimizando (1.15) respecto a r. Se encuentra que ∂u/∂r = 0 en µ ¶1/6 2A12 r0 = σ, Red fcc → 1.09σ. (1.17) A6 √ La constante de la red fcc es a = 2r, el volumen de la celda primitiva es Ω = a3 /4 y la densidad es masa atomica ρ= (1.18) Ω Evaluando la energia potencial en r0 y en r → ∞ se obtiene la energ´ıa de cohesi´on del s´olido ²A2 u0 = − 6 , Red fcc → −8.6². (1.19) 2A12 El modulo de volumen (magnitud inversa de la compresibilidad) B = −V (∂P/∂V )T puede calcularse en terminos de los parametros de Lennard-Jones. A temperatura 0 K, P = −dU/dV = −du/dΩ, luego ∂2u B = Ω 2. (1.20) ∂Ω Procediendo de forma directa se obtiene µ ¶5/2 75² 4² A6 , Red fcc → 3 . B = 3 A12 (1.21) σ A12 σ La estructura de los cristales que existe, es aquella que maximiza la energia de cohesion. Se deja como ejercicio calcular la energia de cohesion las estructuras bcc y sc, y ver que dan menores energias de cohesion que la fcc.

´ ´ 1.2. FUERZAS INTERATOMICAS Y TIPOS DE SOLIDOS

1.2.2

19

Cristales i´ onicos

El ejemplo clasico de estos es el NaCl. El rasgo caracteristico es que estan compuestos por un metal y un nometal situados en las columnas extremas a ambos lados de de la tabla periodica. El elemento metalico tiene una baja energia de ionizacion, mientras que el elemento no metalico tiene alta energia de afinidad, de modo que el sistema pierde energia cuando se trasnfiere un electron del metal hacia el nometal. De este modo, los atomos quedan cargados y entre ellos se establecen fuerzas de tipo Coulomb, mucho m´as intensas que en los s´olidos de gases nobles. Qi Qj φCoulomb (rij ) = (1.22) rij Las fuerzas de atracci´on y repulsi´on Coulombian tienen un efecto neto de hacer colapsar los atomos negativos hacia los positivos, y deben ser contrarestadas por una fuerza repulsiva (de origen cuantico). El modelo mas simple es asociarle un potencial de esfera rigida, o sea ½ 0, si rij > Ri + Rj , (1.23) φrepulsion (rij ) = ∞, si rij < Ri + Rj , donde Ri son los “radios atomicos”. Estos radios han sido determinados mediante un procedimiento empirico de forma tal que las distancias de enlace en un gran numero de compuestos ionicos sea aproximada (por el m´etodo de m´ınimos cuadrados) por la suma de los radios asignados a los atomos en contacto5 . Para simular propiedades dinamicas y termodinamicas a temperaturas no nulas. se usan potenciales suves. Entre los mas utilizados se encuentran el potencial de Buckinham φBuckingham (r) = A exp(−r/ρ) −

C . r6

(1.24)

La energ´ıa de cohesi´on de los cristales ionicos se calcula mediante el mismo m´etodo que para los cristales de gases nobles. Un problema aparece en el c´alculo de la energ´ı a de Coulomb. La serie que define esta energ´ıa es condicionalmente convergente, debido a que la potencia 1/r decrese muy lentamente con r. Las series condicionalmente convergentes pueden ser reordenadas de forma que su suma de igual a cualquier numero real (esto es un teorema del analisis matematico) y ademas su convergencia suele ser lenta. Por lo tanto, obtener de una forma fisicamente correcta y eficiente requiere t´ecnicas especiales, como el m´etodo de Ewald6 .

1.2.3

Cristales covalentes

El paradigma de los cristales covalentes es el diamante. Los cristales compuestos por elementos de la grupo IV y de las columnas vecinas pertenecen a este grupo. Si se examina la estructura de diamante, se observa que el numero de coordinacion es 4 y los cuatro enlaces se disponen formando angulos de 109 grados entre si. Estas estructuras dejan mucho espacio vacio si se consideran los atomos como esferas rigidas, por lo que se llaman estructuras 5 R. D. Shannon, Revised effective ionic radii and systematic studies of interatomic distances in halides and chalcogenides, Acta Cryst. A 32, 751 (1976). 6 P. P. Ewald, Ann. Physik 64, 253 (1921).

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

20

abiertas. Los potenciales de interaccion en los cristales covalentes no pueden ser aproximados mediante ideas de la fisica clasica. La situacion fisica es la misma que en los enlaces covalentes que determinan la estructura de las moleculas en fase gaseosa. Una vez obtenidos de forma numerica los potenciales de interaccion, estos se pueden aproximar por una variedad de formas analiticas. Entre ellos mencionemos el potencial de Morse i h¡ ¢2 (1.25) φM orse = D 1 − e−a(r−r0 ) − 1 . Este potencial es sugerido por la fisica molecular. Para una molecula diatomica r0 es la energia de enlace y D es la energ´ıa de enlace. Para describir los angulos de enlace covelente es necesario tener encuentra potenciales que involucran 3 o mas ´atomos, por ejemplo 1 φ3 (~r1 , ~r2 , ~r3 ) = k2 (θ213 − θ0 )2 exp(−r12 /ρ) exp(−r13 /ρ), 2

(1.26)

donde θ213 es el angulo de enlace covalente centrado en el ´atomo 1. k2 , ρ son par´ametros emp´ıricos que se obtienen mediante ajustes por el m´etodo de m´ınimos cuadrados.

1.2.4

Metales

Al igual que los cristales covalentes, los metales no pueden ser descritos por potential de pares. las propiedades din´amicas resultan mal predichas, adem´as no reproducen la the la llamada discrepancia de Cauchy para las constantes elasticas (C11 6= C44 ). El solo uso d potenciales de pares tambien conduce a estimados incorrectos de las energias de formacion de vacancias, cuyos valores dan muy proximos a las energias de cohesion, mientra que los experimentos indican que deben ser aproximadamente 1/3 de estas. Existen varios tipos de estos potential, todos de muchos cuerpos. Estos potenciales han sido desarrollados para ajustar los calores de constantes de red, energias de cohesion, y constantes elasticas. Como ejemplo, se da la forma de los potenciales de Sutton-Chen · XX X√ ¸ 1 E=² V (rij ) − c ρi , (1.27) 2 i j6=i i donde

µ V (rij ) =

a rij

¶n

X µ a ¶m y ρi = . r ij j6=i

(1.28)

Aqui, rij es la distancia entre los atomos i y j, ² es un parametro con dimensiones de energia, a es un parametro con dimensiones de longitud que normalmentes la constante de la red, c es adimensional, mientras que n y m are positive integers with n > m. Para un cluster diatomico, la distancia de enlace es is dada por ³ m ´1/k m c , k= − n. (1.29) rmin = a n 2 Tambi´en existen cristales que presentan caracter´ısticas mixtas. Un ejemplo cl´asico es el grafito (figura 1.10), que presenta enlaces covalentes muy fuertes en cirtos planos, mientras que la ligaz´on entre los planos se efectua mediante fuerzas del tipo van der Waals.

1.3. ELASTICIDAD

21

Figure 1.10: Estructura del grafito. Los ´atomos de un mismo plano se cohesionan mediante ´ enlaces covalentes. Atomos de distintos planos interact´ uan d´ebilmente por fuerzas de tipo van der Waals. Al escribir con un l´apiz de grafito, se exfolian planos enteros.

1.3

Elasticidad

~ 0 ) de sus posiciones de equlibrio. Consideremos los iones de una red cristalina desplazados ~u(R Si los deslazamientos son peque˜ nos, la energ´ıa potencial del cristal se puede expandir en serie de Taylor hasta los t´erminos de segundo orden U arm = U0 +

1 X ~ µν (R, ~ R ~ 0 )uν (R ~ 0) , uµ (R)D 2 0

(1.30)

∂ 2U ~ ~0 ∂uµ (R)∂u ν (R )

(1.31)

~ R ~ ,µν R,

donde ~ R ~ 0 ) = Dµν (R ~ −R ~ 0) = Dµν (R,

~ −R ~ 0 ) se debe a la simetr´ıa de traslaci´on. Otras propiedades muy La dependencia de D(R importantes que son ~ −R ~ 0 ) = Dνµ (R ~ 0 − R) ~ Dµν (R ~ = Dµν (−R) ~ Dµν (R) X ~ = 0. Dµν (R) ~ R

(1.32) (1.33) (1.34)

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

22

Si las interacciones interat´omicas que se describen para un potencial de pares (ver AshcroftMermin, ec. (22.2-22.11)): U=

1 X ~ ~0 NX ~ φ(R − R ) = φ(R) , 2 2 0 ~ R ~ R,

entonces

(1.35)

~ =0 R6



 ~ −R ~ 0) = δ ~ ~ 0  Dµν (R R,R

X

~ −R ~ 0) , ~ −R ~ 00 ) − φµν (R φµν (R

(1.36)

~ 00 R

donde φµν =

∂ 2 φ(r) , ∂rµ ∂rν

(1.37)

ver demostraci´on en Ashcroft-Mermin cap´ıtulo 22. No obstante, notemos que la ec. (1.30) es util aunque haya potenciales mas generales que los potenciales de pares. Usando las ~ −R ~ 0 ) se puede llevar la propiedades de simetria de la matriz de constantes de fuerzas Dµν (R energia del cristal a la forma U arm − U0 = −

1 X ~ 0 ) − uµ (R)}D ~ ~ ~0 ~0 ~ {uµ (R µν (R, R ){uν (R ) − uν (R)}. 4 0

(1.38)

~ R ~ ,µν R,

A partir de ahora considerarmos la energia del cristal relativa a la energia de equilibrio U0 . En una deformaci´on macroscopica que ocurre cuando un s´olido es sometido a deformaci´on, ~ var´ıan suavemente de una celda a la celda vecina. Entonces podelos desplazamientos uν (R) ~ cuando ~r es un vector de la mos considerar una funci´on continua ~u(~r) que es igual a ~u(R) 0 ~ −R ~ ) se puede hacer la aproximaci´on red de Bravais. Si ~u(~r) var´ıa poco en el rango de D(R ~ 0 ) = ~u(R) ~ + (R ~ 0 − R) ~ · ∇~u| ~ . ~u(R ~ r =R Sustituyendo en la ec. (1.38) si obtiene µ ¶µ ¶ 1 X ∂ ∂ arm ~ ~ U = uµ (R) uν (R) Eσµτ ν , 2 ∂xσ ∂xτ

(1.39)

(1.40)

~ R,µ,ν,σ,τ

donde Eσµτ ν = −

1X ~ τ. Rσ Dµν (R)R 2

(1.41)

~ R

Como las funciones ~u(~r) varian lentamente, se puede escribir (1.40) como una integral ¶µ ¶ µ Z ∂ ∂ 1 X 3 arm uµ (~r) uν (~r) E¯σµτ ν , (1.42) d ~r U = 2 µ,ν,σ,τ ∂xσ ∂xτ donde E¯σµτ ν = Eσµτ ν /Ω, siendo Ω el volumen de la celda primitiva. La ecuaci´on (1.42) es el punto de partida de la teor´ıa macrosc´opica de la elasticidad. El conjunto de 34 = 81 magnitudes E¯σµτ ν es una propiedad de cada material y forma un tensor, transform´andose

1.3. ELASTICIDAD

23

como tal ante rotaciones de los ejes de coordenadas. En el ´algebra de tensores se utiliza el P convenio, debido a Albert Einstein, de omitir los simbolos de sumatioria , entendiendose la suma cada vez que en una formula aparecen dos indices repetidos. As´ı, la ecuaci´on (1.42) se escribe µ ¶µ ¶ Z 1 ∂ ∂ 3 arm d ~r uµ (~r) uν (~r) E¯σµτ ν . U = 2 ∂xσ ∂xτ En lo adelante utilizaremos el arriba mencionado convenio de suma.

1.3.1

Simetr´ıas

De la definici´on (1.41) podemos notar que Eσµτ ν no cambia si se intercambian µ ↔ ν y o τ ↔ σ. Por tanto, es suficiente especificar Eσµτ ν para los siguientes valores de los pares µν y στ xx, yy, zz, yz, zx, xy. (1.43) Esto indica que de las 34 = 81 componentes del tensor Eσµτ ν solo hay 6×6 = 36 componentes independientes. Este n´ umero se reduce mas, si se considera que ante una rotaci´on r´ıgida del cristal la ener´ıa del cristal no se afecta. En una rotaci´on infinitesimal de angulo dω alrededor ˆ , todos los vectores de la red sufren la transformaci´on de un eje de direcci´on ~n ~ −→ R ~ + ~u(R), ~ R

ˆ × R. ~ = δω~n ~ ~u(R)

(1.44)

~ de la ecuacion anterior en (1.40) y exigiendo la que U arm = 0 para δω Sustituyendo ~u(R) arbitrario, se encuentra que U arm s´olo depende de las combinaciones sim´etricas µ ¶ 1 ∂ ∂ εσµ = uµ + uσ . (1.45) 2 ∂xσ ∂xµ

Sistema Triclinico Monoclinico Ortorrombico Tetragonal Romboedrico Hexagonal C´ ubico Amorfo

Table 1.3: Nmero ´ de constantes elasticas independientes. Grupo puntual Constantes elasticas todos 21 todos 13 todos 9 C4 ,C4h ,4 7 C4v , D4 , D4h , D2d 6 C3 , S6 7 C3v , D3 , D3d 6 todos 5 todos 3 2

El tensor sim´etrico εσµ se denomina tensor de deformaci´ on. Consecuentemente, se puede rescribir (1.42) como Z 1 arm d3~rεσµ cσµτ ν ετ ν , (1.46) U = 2

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

24 donde cσµτ ν = −

1 X [Rσ Dµν Rτ + Rµ Dσν Rτ + Rσ Dµτ Rν + Rµ Dστ Rν ] . 8Ω

(1.47)

~ R

De (1.47) se deduce que cσµτ ν es invariante ante las permutaciones σµ ↔ τ ν, σ ↔ µ y τ ↔ ν. Como resultado, el n´ umero de componentes independientes de reduce a 21. Las propiedades anteriores son generales y v´alidas para cualquier sistema cristalino. Se puede reducir mas el numero de constantes elasticas independientes en dependencia del grupo puntual de simetria del cristal y el tipo de red de Bravais. La tabla 1.3 resume el n´ umero de constantes el´asticas para todos los sistemas cristalinos y para los amorfos. Por ejemplo, en el caso c´ ubico, las u ´nicas tres componentes independientes son C11 = cxxxx = cyyyy = czzzz C12 = cxxyy = cyyzz = czzxx C44 = cxyxy = cyzyz = czxzx .

(1.48) (1.49) (1.50)

Todas las demas componentes, en las cuales x, y o z aparece un numero impar de veces, son 0. En las ecuaciones anteriores se usa el convenio xx ≡ 1,

yy ≡ 2,

zz ≡ 3,

yz ≡ 4,

zx ≡ 5,

xy ≡ 6.

(1.51)

Convencionalmente, en teoria de la elasticidad se utiliza una notaci´on ligeramente modificada. El campo de desplazamiento se describe por la magnitud llamada deformaci´ on, relacionada con el tensor de deformaci´on seg´ un eµν = εµν , si µ = ν = 2εµν , si µ 6= ν,

(1.52) (1.53)

cuya notaci´on se simplifica a ei = eµν de acuerdo a (1.51). En lugar de (1.46) se escribe 6 Z 1X U= d3~rei Cij ej , 2 i,j=1

(1.54)

donde Cij = cσµτ ν , acorde a (1.51). Las cantidades Cij forman una matriz de dimension 6 × 6 (no es un tensor) y se denominan modulos el´ asticos (elastic moduli o stiffness constants). Los elementos de la matriz S que es inversa a C se denominan constantes el´ asticas (elastic constants o elastic compliance constants). Una aplicaci´on de la teor´ıa de la elasticidad es la ecuaci´on de las ondas el´asticas. La energ´ıa cin´etica asociada a un campo de deformaci´on es Z 1 ˙ ρ~u(~r, t)2 d3~r, (1.55) T = 2 donde ρ es la densidad. El Lagrangiano del medio es ¸ Z · 1 ˙2 1 L=T −V = ρ~u − εσµ cσµτ ν ετ ν d3~r. 2 2

(1.56)

1.3. ELASTICIDAD

25

El principio variacional de Hamilton Z Z δ Ldt = δ L d3~rdt = 0 conduce a las ecuaciones de Lagrange µ ¶ µ ¶ ∂ ∂ ∂L ∂L ∂L − − =0 ∂uµ ∂t ∂ u˙ µ ∂xν ∂uµ,ν

µ ∂uµ u˙ µ = , ∂t

(1.57)

uµ,ν

∂uµ = ∂xν

¶ .

(1.58)

Haciendo aproximadamente unas 3 p´aginas de algebra las ecuaciones de Lagrange se reducen a la forma ∂uτ ∂ ρ¨ uµ = cµσντ = (cµσντ εντ ) . (1.59) ∂xσ ∂xν ∂xσ Tarea. Demuestre la ecuaci´on (1.59). El tensor σµσ = cµσντ εντ (1.60) se denomina tensor de esfuerzos y la ecuacion (1.60) es la Ley de Hooke. ρ¨ uµ es la fuerza por unidad de volumen y segun la ecuaci´on (1.59), es igual a la divergencia del tensor de esfuerzos. Esto es plenamente consistente con el teorema de Gauss para tensores Z Z I ∂σµσ fµ dV = dV = σµσ dSσ . (1.61) ∂xσ La ecuaci´on anterior expresa el hecho fisico de que las fuerzas internas se anulan y la fuerza total es igual a la suma de las fuerzas aplicadas en la superficie. El elemento σµσ dSσ

(1.62)

~ (area dS es la componente µ de la fuerza que act´ ua sobre el elemento de superficie dS y direccion perpendicular dada por el vector). Tomando elementos de superficie en los planos xy,yz,zx, encontramos que la componente σµσ es la componente µ de la fuerza que actua sobre la unidad de area perpendicular al eje xσ . Por ejemplo, para un area dS paralela al plano yz, ~ debe ser perpendicular al plano yz: dS1 = dS, ds2 = dS3 = 0. La tenemos que el vector dS fuerza total ejercida en el area dS es la suma de la fuerza normal y la fuerza tangencial: σ11 es la fuerza normal al plano (como una presion hidrostatica), σ21 y σ31 son las componentes de la fuerza tangencial (como la fuerza de roce). Con esto se describen todas las formas posibles de fuerzas actuantes en una superficie. La forma σµσ dSσ (sumado sobre σ seg´ un convenio de Einstein) es la expresion general cuando el plano no es paralelo a uno de los ejes coordenados (o cuando el sistema de referencia no tiene un eje paralelo al plano). Debido a las propiedades de simetr´ıa del los tensores εµν y cµσντ , el tensor de esfuerzos es un tensor sim´etrico y por tanto tiene 6 magnitudes independientes. Estas se reunen en una matriz columna de dimension 6, ti = σµν , siguiendo el convenio dado por la f´ormula (1.51), (note que (1.52) y (1.53) no se aplican a la tensi´on). A este vector columna se le llama tensi´ on (note que no es un tensor, el tensor es σµσ ). Con esta notaci´on, la Ley de Hooke (1.60) se rescribe 6 X ti = Cij ej . (1.63) j=1

N´otese los las matrices ti , Cij y ej no son tensores, pues no transforman como tales ante transformaciones cartesianas. Por eso se ha escrito explicitamente el signo de suma.

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

26

1.3.2

Deformaciones homog´ eneas

Consideremos un ortoedro microsc´opico de lados x1 , x2 , x3 en un cuerpo sin deformar. En un sistema de coordenadas adecuado con origen en un vertice del ortoedro, los demas vertices se definen en funcion de los vectores xˆ = (x1 , 0, 0), yˆ = (0, y2 , 0) y zˆ = (0, 0, z3 ). Consideremos ahora una deformacion homogenea descrita por el tensor εν . Los vectores cambian a xˆ0 = (x1 + u1 , u2 , u3 ), yˆ0 = (v1 , y2 + v2 , v3 ), zˆ0 = (w1 , w2 , z3 + w3 ).

(1.64) (1.65) (1.66)

Si la deformaci´on es peque˜ n a entonces u1 , u2 , u3 son proporcionales a x1 , v1 , v2 , v3 son proporcionales a y2 y w1 , w2 , w3 a z3 . Los factores de proporcionalidad son las componentes del tensor de deformaciones ∂u2 ∂u3 ∂u1 x1 = ε11 x1 , u2 = x1 = ε˜12 x1 u3 = x1 = ε˜13 x1 ∂x1 ∂x1 ∂x1 ∂v2 ∂v1 = y2 = ε22 y2 , v1 = y2 = ε˜21 y2 , etc, ∂y2 ∂y2 ∂w3 = z3 = ε33 z3 , etc. ∂z3

u1 =

(1.67)

v2

(1.68)

w3

El volumen inicial del ortoedro es V = x1 y2 z3 . Al ser deformado es ¯ ¯ ¯ x 1 + u1 ¯ u u 2 3 ¯ ¯ 0 ¯ = x1 x2 x3 + u1 y2 z3 + v2 x1 z3 + w3 x1 y2 + ... v1 y2 + v2 v3 V = ¯¯ ¯ ¯ w1 w2 z3 + w3 ¯

(1.69)

(1.70)

Los t´erminos de orden cuadr´atico o c´ ubico en los productos uvw son despreciables si las componentes del tensor de deformaci´on son peque˜ nas y µ ¶ u1 v 2 w 3 ∆V = + + = ε11 + ε22 + ε33 = Trε = e1 + e2 + e3 . (1.71) V x1 y2 z3 La ecuaci´on anterior dice que el cambio relativo de volumen es igual a la traza del tensor de deformaci´on, o equivalentemente, la suma de las tres primeras componentes de la deformaci´ on. Adem´as, las componentes ε11 , ε22 , ε33 indican la contracci´on o dilataci´on relativa a lo largo de cada eje coordenado. Los t´erminos no diagonales tambi´en tienen interpretaci´on. Consideremos el ´angulo entre los vectores deformados, por ejemplo, cos(ˆ x0 , yˆ0 ) =

x1 v1 + y2 u2 xˆ0 · yˆ0 ' = ε˜12 + ε˜21 = 2ε12 . 0 0 |ˆ x ||ˆ y| x1 y 2

(1.72)

En la secuencia anterior se despreciaron los productos de dos o mas componentes de las deformaciones u, v, w. As´ı escribimos εxy =

1 cos(ˆ x0 , yˆ0 ), 2

εyz =

1 cos(ˆ y 0 , zˆ0 ), 2

εzx =

1 cos(ˆ z 0 , xˆ0 ). 2

(1.73)

1.3. ELASTICIDAD

27

De la ecuaci´on anterior se observa que si los ´angulos del ortoedro se conservan en 90◦ , el tensor de deformaci´on (referido a esos mismos ejes) es diagonal. Una deformaci´on en la cual el cambio de volumen es nulo, se llama deformaci´ on de corte. En este caso el tensor de deformaci´on tiene la forma general 1 = τµν − (Trτ )δµν . εcorte µν 3

(1.74)

Una deformac´on que cambia el volumen pero no la forma del cuerpo, se llama deformaci´ on hidrost´ atica y su tensor de deformaci´on tiene la forma εhidro = constante × δµν . µν

(1.75)

Cualquier deformaci´on se puede representar como la suma de una deformaci´on hidrostatica y una de corte, gracias a la identidad µ ¶ 1 1 + εhidro . (1.76) εµν = εµν − (Trε)δµν + (Trε)δµν = εcorte ν ν 3 3 Examinemos como cambia la energia de la red ante una deformacion hidrostatica en el caso particular de un cristal c´ ubico. La deformaci´on es e1 = e2 = e3 = ε, e4 = e5 = e6 = 0. Usando la ec. (1.54) para una deformaci´on homogenea obtenemos δU 1X = ei Cij ej V 2 i,j    C11 C12 C12 0 0 0 ε  C12 C11 C12 0   0 0  ε     ε  ¡ ¢  C12 C12 C11 0 0 0   ε ε ε 0 0 0  = (1.77)  0   0 0 C44 0 0    0   0 0 0 0 C44 0   0  0 0 0 0 0 C44 0 µ ¶2 ∆V 1 1 (C11 + 2C12 )(3ε)2 = (C11 + 2C12 ) . (1.78) = 6 6 V Comparando con la definici´on del modulo de volumen (1.20) se obtiene que 1 B = (C11 + 2C12 ). 3

(1.79)

Un caso importante en ingenieria civil cuando sobre un cuerpo actua un esfuerzo unidireccional. Considerese una viga de acero sometido a una presi´on vertical (eje Z). Debido a este esfuerzo la viga se deforma, expandiendose un poco en las direcciones laterales. Apliquemos la ley de Hooke. La tensi´on es t3 = σzz = P , t1 = t2 = t4 = t5 = t6 = 0.      e1 C11 C12 C12 0 0 0 0    0   C12 C11 C12 0 0 0    e2        P   C12 C12 C11 0 0 0   e3    = (1.80)   e4  .  0   0 0 0 C 0 0 44       0   0 0 0 0 C44 0   e5  e6 0 0 0 0 0 C44 0

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

28

Resolviendo la ecuaci´on anterior se obtiene c12 P , (C11 − C12 )(C11 + 2C12 ) (C11 + C12 )P = , (C11 − C12 )(C11 + 2C12 ) = 0.

e1 = e2 = − e3 e4 = e5 = e6

(1.81) (1.82) (1.83)

La relaci´on (1.82) sirve de base para definir el M´ odulo de Young E=

P (C11 − C12 )(C11 + 2C12 ) = . e3 (C11 + C12 )

(1.84)

Las relaciones (1.81) y (1.82) justifican la definici´on del Coeficiente de Poisson σ=−

e1 εxx C12 =− = . e3 εzz C11 + C12

(1.85)

El m´odulo de Young E y el coeficiente de Poisson σ son magnitudes c´omodas para los c´alculos de resistencia de materiales habituales el ingenier´ıa. Adem´as, en medios isotropos solamente hay dos constantes el´asticas independientes Cij 7 , las cuales se expresan en funci´on de E y σ. Las formulas pertinentes pueden verse en Landau y Lifshitz, Teoria de la Elasticidad. En particular, el m´odulo de volumen es en un medio c´ ubico o is´otropo es igual a B=

1.3.3

E . 3(1 − σ)

(1.86)

Ondas el´ asticas

Consideremos las soluciones de la ecuaci´on (1.59) que dependen del tiempo. En un medio infinito y homog´eneo se puede proponer la soluci´on en forma de onda plana ~u(~r, t) = ~² exp[i(~k · ~r − ωt)], se obtiene la ecuaci´on de autovalores ρω 2 ²µ =

à X X τ

(1.87)

! cµσντ kσ kν

²τ .

(1.88)

στ

Para mayor claridad hemos mantenido el signo de suma en la ecuaci´on anterior. Esta es la ecuaci´on que describe la propagaci´on del sonido en un cristal anisotropo. Consideremos la ecuacion anterior para el caso de un cristal cubico. Utilizando las relaciones (1.50) se obtienen las siguientes ecuaciones para el vector ~²      ²1 (C12 + C44 )k1 k3 C11 k12 + C44 (k22 + k32 ) (C12 + C44 )k1 k2 ²1   ²2  C11 k22 + C44 (k12 + k32 ) (C12 + C44 )k2 k3 ρω 2  ²2  =  (C12 + C44 )k1 k2 2 2 2 ²3 (C12 + C44 )k1 k3 (C12 + C44 )k2 k3 C11 k3 + C44 (k1 + k2 ) ²3 (1.89) 7

En la pr´oxima secci´on obtendremos una f´ormula para C44 en funci´on de C11 y C12 en un medio is´otropo.

1.3. ELASTICIDAD

29

La ecuaci´on anterior es un problema de autovalores para cada valor especificado del vector de onda ~k = (k1 , k2 , k3 ). En algunos casos es simple obtener soluciones analiticas. a) ~k = (k1 , 0, 0). El sistema de ecuaciones queda completamente desacoplado (C11 k12 − ρω 2 )²1 = 0, (C44 k12 − ρω 2 )²2 = 0, (C44 k12 − ρω 2 )²3 = 0.

(1.90) (1.91) (1.92)

Las soluciones son

p p a.1) ²1 6= 0, ²2 = ²3 = 0 (modo longitudinal), ω = C11 /ρk1 , vL = C11 /ρ. p p a.2) ²2 6= 0, ²1 = ²3 = 0 (modo transversal), ω = C44 /ρk1 , vT = C44 /ρ. p p a.3) ²3 6= 0, ²1 = ²2 = 0 (modo transversal), ω = C44 /ρk1 , vT = C44 /ρ. Las ondas el´asticas son las ondas sonoras, entendiendo que la frecuencia ω/2π est´e en el rango audible 5 − 20000 Hz. Midiendo las velocidades de las ondas transversales y longitudinales a lo largo de uno de los ejes principales del cristal, se determinan las constantes el´asticas C11 y C44 .√ b) ~k = (k/sqrt2, k/ 2, 0). En este caso el sistema no se desacopla totalmente, quedando      2 2 ²1 (C11 + C44 ) k2 (C12 + C44 ) k2 0 ²1 2 2   ²2  . ρω 2  ²2  =  (C12 + C44 ) k2 (C11 + C44 ) k2 0 (1.93) 2 ²3 ²3 0 0 C44 k Las soluciones son

p

(C11 + C12 + 2C44 )/2ρ k = vL k. p b.2) ²2 = −²1 6= 0, ²3 = 0 (modo transversal), ω = (C11 − C12 )/2ρ k = vT k k. p p b.3) ²3 6= 0, ²1 = ²2 = 0 (modo transversal), ω = C44 /ρ k, vT ⊥ = C44 /ρ.

b.1) ²1 = ²2 6= 0, ²3 = 0 (modo longitudinal), ω =

Las relciones anteriores permiten determinar la otra constante el´astica. N´otese que la velocidad de la onda depende de la direcci´on de propagaci´on. Esto es una manifestaci´on de la anisotrop´ıa de los cristales. En un medio el´astico is´otropo la velocidad de propagaci´on no puede depender de la direcci´on de propagaci´on. Si uno iguala las velocidades transversales de la soluci´on b) entre si y con las de la solucion (a), se llega a la relaci´on 2C44 = C11 − C12

en un medio is´otropo.

(1.94)

Si se igualan las velocidades longitudinales de las soluciones (a) y (b), se obtiene el mismo resultado. Los medios el´asticamente is´otropos son materiales amorfos, policristales (por ejemplo, el acero), o algunos cristales muy especiales. Las constantes el´asticas de un policristal no son exactamente el promedio de las constantes del monocristales, como podria pensarse ingenuamente. Para determinarlas te´oricamente ser´ıa necesario resolver el sistema de ecuaciones de la Ley de Hooke para el conjunto de todos los granos que forman el policristal.8 8 Si alguien encuentra una referencia en que se reporte haber hecho este c´alculo, agradecere que me la muestre.

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

30

1.4

Modelos discretos de las vibraciones de la red cristalina

En esta secci´on examinaremos los l´ımites de la teor´ıa de la elasticidad, debido al hecho de que la distribuci´on de masa no es realmente continua. Para simplificar el tratamiento matem´atico sin perder la esencia del fen´omeno, consideremos un modelo de movimiento unidimensional, que a pesar de su simplificaci´on, tiene utilidad pr´actica.

1.4.1

Modos normales de una cadena lineal monoat´ omica

Figure 1.11: Deformaciones de una red cristalina unidimensional de par´ametro a. Consideremos que los iones vecinos interact´ uan el´asticamente 1 X U arm = K [u(na) − u[(n + 1)a]]2 . 2 n

(1.95)

Ecuaci´on de movimiento 1 ∂U arm = − K [2(u(na) − u((n + 1)a)) − 2(u((n − 1)a) − u(na))] ∂u(na) 2 = −K[2u(na) − u((n + 1)a) − u((n − 1)a)] .

M u¨(na) = −

(1.96) (1.97)

Si la cadena tiene longitud N (i.e. iones en a, 2a, · · · , N a) las ecuaciones din´amicas se completan con las condiciones de borde. Postulado: Si N → ∞, lo que ocurre al interior no depende de las condiciones de borde. La forma m´as conveniente es u(0) = u(N a), conocida como condici´on de Born-von Karman o condici´on de borde peri´odica. Consideremos soluciones de la forma u(na, t) = ²ei(kna−ωt) ,

(1.98)

con la condici´on de borde que eikN a = 1. La condici´on anterior implica kN a = 2π` , 0

` ∈ Z, 2π

k=

2π l . a N

(1.99)

N´otese que k 0 = k + 2π , eik na = ei a ·na eikna = eikna . Entonces k 0 y k dan exactamente a la misma soluci´on. Luego, las soluciones independientes pueden concentrase en cualquier

1.4. MODELOS DISCRETOS DE LAS VIBRACIONES DE LA RED CRISTALINA

31

£ ¤ intervalo de longitud 2π/a. Es usual considerar el intervalo k ∈ − πa , πa . Tomando en cuanta la ec. 1.99, expresaremos k en la forma k=

2π ` , a N

`=−

N N + 1, · · · , . 2 2

(1.100)

Reemplazando (1.98) en la ecuaci´on de movimiento se obtiene −M ω 2 ei(kna−ωt) = −K[2 − e−ika − e−ika ]ei(kna−ωt) = −2K(1 − cos ka)ei(kna−ωt) De donde se obtiene la relaci´on r r 2K(1 − cos ka) K ω = ω(k) = =2 M M

¯ ¯ ¯ ka ¯ ka π π ¯sin ¯ , ∈ [− , ] . ¯ ¯ 2 2 2 2

(1.101)

(1.102)

π/a

-π/a

Figure 1.12: Ley de dispersi´on de los modos normales de vibraci´on en una cadena lineal monoat´omica. Las soluciones reales independientes son las partes reales o imaginarias de ²ei(kna−ωt) uk (na, t) = Ak cos (kna − ωt) + Bk sin (kna − ωt),

(1.103)

p donde ω = 2K/M | sin(ka/2)|. Tenemos N valores de k, luego, hay 2N constantes Ak y Bk . Un movimiento arbitrario se determina totalmente por N posiciones iniciales y N velocidades iniciales. X U general (na, t) = Ak cos (kna − ωt) + Bk sin (kna − ωt) . (1.104) k

Evaluando en t = 0 u(na, 0) =

X

Ak cos kna + Bk sin kna

k

u(na, ˙ 0) =

X

{+ωAK sin kna − ωBk cos kna} .

k

Las condiciones iniciales determinan las constantes Ak , Bk .

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

32 Vibraciones de onda larga

¿Qu´e significa exactamente que k sea peque˜ no o que la longitud de onda sea grande? El criterio es comparar la longitud de onda con la longitud caracter´ıstica a. Esto es equivalente a la condici´on ka ¿1 2 Ã r ! K |k| , (1.105) ω(k) = a M q K la velocidad del sonido. Esto es caracter´ıstico de un medio homog´eneo, pues es siendo a M lo que se obtiene de las ecuaciones macroscopicas de la teor´ıa de la elasticidad. Cuando ka es del orden de 1, se deja de cumplir la linealidad en la ley de dispersi´on. Principio f´ısico: La no linealidad aparece cuando λ ∼ a, siendo a la longitud caracter´ıstica de la inhomogeneidad. Nota: hay soluciones con ω < 0. Si se cambia ω por −ω y k por −k se obtiene la misma soluci´on para ei(kn−ωt) .

1.4.2

Redes 1-D con bases

Consideremos una red unidimensional con 2 iones por celda primitiva

Figure 1.13: Red unidimensional con 2 iones por celda primitiva.

U arm = U0 +

KX G [u1 (na) − u2 (na)]2 + [u2 (na) − u1 ((n + 1)a)]2 2 n 2

(1.106)

Ecuaciones de movimiento: ∂U arm = −K[u1 (na) − u2 (na)] − G[u1 (na) − u2 ((n − 1)a)] ∂u1 (na) ∂U arm M u¨2 (na) = − = −K[u2 (na) − u1 (na)] − G[u2 (na) − u1 ((n + 1)a)] ∂u2 (na)

M u¨1 (na) = −

(1.107) (1.108)

La soluci´on propuesta es u1 (na) = ²1 ei(kna−ωt) ,

u2 (na) = ²2 ei(kna−ωt)

(1.109)

1.4. MODELOS DISCRETOS DE LAS VIBRACIONES DE LA RED CRISTALINA

33

Sustituyendo en las ecuaciones de movimiento (1.107) y (1.108) se obtiene [M ω 2 − (K + G)]²1 + (K + Ge−ika )²2 = 0 (K + Geika )²1 + [M ω 2 − (K + G)]²2 = 0 Este sistema de ecuaciones tiene soluci´on no trivial si ¯ ¯ −ika ¯ ¯M ω 2 − (K + G) K + Ge ¯=0 ¯ −ika 2 ¯ K + Ge M ω − (K + G)¯

(1.110)

Siendo la soluci´on ω2 =

K +G 1√ 2 ± K + G2 + 2KG cos ka , M M

²2 K + Ge−ika =∓ . ²1 |K + Ge−ika |

(1.111)

Hay cuatro ω. Basta con tener las positivas, pues ei(kna−(−ω(k)t)) = ei(kna+ω(k)t) = ei(kna+ω(−k)t) .

(1.112)

Las soluciones reales son: cos (kna + ω(−k)t) = cos ((−k)na − ω(−k)t) sin (kna + ω(−k)t) = − sin ((−k)na − ω(−k)t) ,

(1.113)

Es decir, tomar −ω(~k) equivale a la onda con −~k, ω(−~k). En resumen: ´ Para todo k hay dos ω(~k) y dos parejas (²1 , ²2 ). La rama inferior se denomina ACUSTICA ´ y la rama superior se denomina OPTICA

Figure 1.14: Relaci´on de dispersi´on de fonones para una red unidimensional con dos iones en la celda primitiva. Notemos que presenta dos ramas: la inferior o ac´ ustica y la superior u ´optica. Tenemos el caso donde ka ¿ 1 de manera que cos ka ∼ 1 −

(ka)2 luego 2

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

34 s ω= r ω=

KG ka , 2M (K + G)

²1 = ²2

2(K + G) − O(ka)2 , 2

ac´ ustica movimiento en fase,

²2 = −²1

´optica movimiento en antifase.

(1.114) (1.115)

Figure 1.15: Esquema del movimiento de los iones en la cadena diatomica lineal para ondas largas (k ∼ 0). a) Rama acustica. b) Rama ´optica. Como muestra la Fig. 1.15, en los modos ac´ usticos de onda larga los ´atomos de la misma celda oscilan en fase, mientras que en el modo ´optico lo hacen en contrafase. Como veremos m´as adelante, esto es un caso particular de una propiedad general que establece que en los modos ´opticos de onda larga, el centro de masas de la celda elemental permanece inmovil. Consideremos ahora el l´ımite opuesto, cuando k = πa , K > G r 2G ω= , ²1 = +²2 ac´ usticas (1.116) M r 2K ω= , ²1 = −²2 ´opticas (1.117) M

Figure 1.16: Esquema del movimiento de los iones en la cadena diatomica lineal para ondas cortas (k ∼ π/a). a) Rama acustica. b) Rama ´optica. N´otese que en el caso k = π/a, se cumple que eikna = einπ y la fase cambia de sitio en sitio. Si ²2 = −²1 , coinciden las fases de los iones consecutivos si pertenecen a distintas celdas. Si el caso es que K = G se cierra el gap. En este caso debe ocurrir que la constante ustica solamente. El pnto de vista con de la red es a2 y en esta red tendremos una rama ac´ constante a corresponde a doblar la zona de Brillouin.

´ DE ESTADO 1.5. ECUACION

35

Figure 1.17: Relaci´on de dispersi´on para una red unidimensional con dos iones en la celda primitiva, y iguales constantes el´asticas.

1.5

Ecuaci´ on de estado

En los cursos elementales de f´ısica se aprende la famosa ecuaci´on de estado del gas ideal pV = N kB T.

(1.118)

Su deducci´on se fundamenta en que las part´ıculas tienen dimensiones despreciables (comparadas con las distancias entre ellas) y que interact´ uan entre s´ı y con las paredes del recipiente solo lo hacen mediante choques instant´aneos. Para un sistema de part´ıculas interactuantes en equilibrio termodin´amico existe en principio una ecuaci´on de estado, si bien mas complicada y dificil de obtener.

1.5.1

Ecuaci´ on de Murnaghan

Consideremos un medio s´olido sometido a una compresi´on hisdrost´atica, cuyo modulo de volumen se define y aproxima por ¶ µ ∂p = B0 + B00 p. (1.119) B = −V ∂V T Comentemos la ecuaci´on anterior. En teoria de elasticidad hemos obtenido expresiones para B en t´erminos de las constantes el´asticas, por ejemplo (1.79) y (1.86). Estas f´ormulas corresponden a B0 cuando la presi´on externa es 0 (estrictamente p ¿ B/B00 , lo cual se cumple a la presi´on atmosf´erica). Ahora, cuando las deformaciones dejan de ser pequeas se hace necersario considerar la dependencia con la presi´on externa (o con el cambio de volumen) y B00 es simplemente el t´ermino de primer orden en la expansi´on de McLaurin. La ecuaci´on diferencial (1.119) puede integrarse a temperatura constante, obteniendo la relaci´on "µ ¶ 0 # B0 B0 V0 p= 0 −1 . (1.120) B0 V Esta es la Ecuaci´on de Estado de Murnaghan.9 Si se puede determinar experimentalmente o teoricamente la dependencia de p con V , un ajuste de minimos cuadrados permite determinar B0 y B00 , que en principio dependen de la temperatura. En las determinaciones teoricas se 9 F. D. Murnaghan, The compressibility of media under extreme pressures, Proc. Nat. Acad. Sci. (PNAS) 30, 244 (1944).

´ CHAPTER 1. ELEMENTOS DE F´ISICA DE SOLIDOS

36

Figure 1.18: Curva P-V para la aleaci´on Ti3 Si0.5 Ge0.5 C2 . Manoun et al, Appl. Phys. Lett. 84, 2799 (2004). calcula directamente la energ´ıa o la energ´ıa libre. Para temperaturas distintas de 0 K la funci´on termodin´amica minimal a V y T constantes es la energ´ıa libre F y p = − (∂F/∂V )T . Integrando la presi´on se obtiene la ecuaci´on de Murnaghan para la energ´ıa libre " # µ ¶B00 −1 B0 V0 V0 V 1 B00 F (V ) = F0 + + − . (1.121) B00 B00 − 1 V V0 B00 − 1 Cualquiera de las dos ecuaciones puede usarse para ajustar los valores de B0 y B00 , ademas de F0 en el caso de (1.121), si bien F0 no tiene un significado absoluto. En simulaciones computacionales a veces el calculo de la energia libre es muy costoso, y se reemplaza por la energ´ıa. Esto es estrictamente valido solamente a temperatura 0 K. El modulo de volumen a presi´on no nula se calcula a partir de (1.119) y (1.120) µ B(V ) = B0

1.5.2

V0 V

¶B00 .

(1.122)

Ecuaci´ on de Birch-Murnaghan

Un refinamiento a la ecuaci´on de Murnaghan fue dado por Birch.10 El resultado se conoce como ecuaci´on de estado de Birch-Murnaghan "µ ¶ "µ ¶ #) µ ¶5/3 # ( 7/3 2/3 V0 V0 3 0 V0 3 − −1 1 + (B0 − 4) , (1.123) p(V ) = B0 2 V V 4 V 10

Birch, Phys. Rev. 71, 809 (1947), J. Geophys. Res. 83, 1257 (1978).

1.6. DENSIDAD DE ESTADOS y

37

" # 7/3 5/3 3 9 V V V E(V ) = E0 − B0 (4 − B00 ) 02 − (14 − 3B00 ) 04/3 + (16 − 3B00 ) 02/3 . 16 V V V

1.6

(1.124)

Densidad de Estados

Para un sistema cualquiera con un espectro de energ´ıas En , la densidad de estados se define X δ(E − En ), (1.125) ρ(E) = n

siendo n el conjunto de numeros que definen un estado del sistema. Si se integra Z E2 X Z E2 X ρ(E)dE = δ(E − En )dE = 1 E1

n

E1

(1.126)

n

E1 Emax el atomo se ioniza. Este campo maximo es del orden de e/r0 . Eval´ uelo para un ´ atomo de H y comp´ arelo con los campos que se miden en el laboratorio. La condici´on de equilibrio ante un campo externo E, para x/r0 ¿ 1 es E = Erest = −

ex ex2 + . 6r02 8r04

(2.6)

63 Despejando x en funci´on del campo externo se obtiene µ 3¶ ¶ µ 6r0 27r05 x= E 2. E+ 2 8 e

(2.7)

El momento dipolar es µ = ex =

27p0 6p0 E + 2 E 2 + ..., E0 E0

p0 = er0 ,

E0 = e/r02 .

(2.8)

La relaci´on anterior se escribe en forma general, para un medio con N/V atomos por volumen P =

Nµ = χ(1) E + χ(2) E 2 + .... V

(2.9)

Incluso en la luz laser se cumple que E ¿ E0 , de modo que χ(2) ¿ χ(1) . En cristales las propiedades no son is´otropas y deben ser expresadas por medio de tensores. Expandiendo hasta tercer orden se escribe (1)

(2)

(3)

Pi = χij Ej + χijk Ej Ek + χijkl Ej Ek El .

(2.10)

Muchos cristales, as´ı como los medios is´otropos, tienen simetr´ıa de inversi´on y para ellos (2) el tensor χijk es nulo. 2 Veamos las implicaciones de la respuesta no lineal. Para eso consideremos las ecuaciones de Maxwell en un medio no lineal sin cargas libres y sin anisotrop´ıa. ~ = 0, ∇·D ~ ~ = − 1 ∂B , ∇×E c ∂t

~ = 0, ∇·B

(2.11)

~ ~ = 1 ∂D ∇×H c ∂t

(2.12)

junto con las relaciones constitutivas del medio ~ = H, ~ B

~ ij Ej = ²(E)Dj . Di = ²(E)

(2.13)

Combinando las ecuaciones de Maxwell se obtiene ~+ ∇×∇×E

~ 1 ∂2E 4π ∂ 2 P~ = − c ∂t2 c2 ∂t2

(2.14)

~ y la primera ecuaci´on de Mawell ~ = ∇(∇ · E) ~ − ∇2 E Usando la identidad ∇ × ∇ × E

obtenemos que

~ = ∇ · (²E) ~ = ²∇ · E ~ +E ~ · (∇²) 0=∇·D

(2.15)

~ = −1E ~ · ∇(ln ²) ' 0. ∇·E ²

(2.16)

2 Recu´erdese que un tensor cartesiano Aijk se transforma ante propiedades de simetria como el producto de coordenadas xi xj xk . La inversi´on cambia xi xj xk por (−xi )(−xj )(−xk ) = −xi xj xk y si es una transformaci´on de simetr´ıa, entonces Aijk = −Aijk .

64

´ CHAPTER 2. ELEMENTOS DE OPTICA NO LINEAL

El caso t´ıpico es que ln ² var´ıa lentamente, y por eso se puede despreciar la divergencia del vector de intensidad de campo el´ectrico. Separando el vector de polarizaci´on en sus partes ~ + 4π P~L = ²L E, ~ podemos llevar la ecuaci´on lineal y no lineal P~ = P~L + P~N L , y utilizando E (2.14) a la forma 2 2~ 2~ ~ − n ∂ E = − 4π ∂ PN L , ∇2 E (2.17) c2 ∂t2 c2 ∂t2 √ donde n = ²L es el ´ındice de refracci´on de la luz. La ec. (2.17) se puede resolver por un m´etodo iterativo, siempre que los efectos no lineales sean correcciones a la soluci´on lineal. Se puede reescribir (2.17) de forma abstracta, como ~ = S(E), ~ L(E)

(2.18) 2

2

n ∂ , c2 ∂t2 2~ ~ = − 4π ∂ PN L S(E) c2 ∂t2 L = ∇2 −

(2.19) (2.20)

~ 0 (que puede ser la soluci´on de la ecuaci´on lineal) y y se resuelven a partir de una semilla E siguiendo un proceso iterativo ilustrado en la ecuaci´on (2.21) n o n o ~ ~ ~ ~ ~ ~ ~ ~ E0 −→ S(E0 ) −→ L(E1 ) = S(E0 ) −→ E1 −→ S(E1 ) −→ L(E2 ) = S(E1 ) −→ ... (2.21) Consideremos el caso de una onda monocromatica de frequencia ω que incide en un medio cuya susceptibilidad de segundo orden χ(2) es no nula. Considerando el medio is´otropo para ~ de modo que simplicar, podemos orientar el eje de coordenadas a lo largo del vector E, trabajaremos solo con ecuaciones escalares E0 = E(ω)eiωt + E(ω)∗ e−iωt .

(2.22)

¢2 ¡ PN L = χ(2) E(ω)eiωt + E(ω)∗ e−iωt = PN L (0) + PN L (2ω)ei2ωt + c.c.,

(2.23) (2.24)

Donde c.c. significa complejo conjugado y PN L (0) = 2χ(2) |E(ω)|2 , PN L (2ω) = χ(2) E(ω)2 .

(2.25) (2.26)

El t´ermino PN L (0) es un campo est´atico generado por un efecto no lineal. El t´ermino PN L (2ω) es el responsable del efecto llamado generaci´ on del segundo arm´ onico. La generaci´on del segundo arm´onico se utiliza en la generaci´on de luz laser. Por ejemplo, el l´aser de rub´ı tiene una longitud de onda de 696 nm (infrarrojo). Al atravesar un cristal de un material denminado KDP, se genera luz de 347 nm, ademas de luz de 696 nm.

´ 2.1. EFECTO ELECTROOPTICO

2.1

65

Efecto electro´ optico

Consideremos un medio sometido a una onda luminosa linealmente polarizada y a un campo electrico constante paralelo al campo de la onda luminosa. Einc = E(0) + E(ω)eiωt + c.c., 2 PN L = χ(2) Einc .

(2.27) (2.28)

Efectuando la operaci´on de elevar al cuadrado obtenemos PN L = PN L (0) + PN L (ω)eiωt + PN L (2ω)ei2ωt + c.c., ¡ ¢ PN L (0) = χ(2) E(0)2 + 2|E(ω)2 | , PN L (ω) = 2χ(2) E(0)E(ω), PN L (2ω) = χ(2) E(ω)2 .

(2.29) (2.30) (2.31) (2.32)

La polarizaci´on lineal se ve afectada por PN L (ω) £ ¤ P (ω) = P (1) (ω) + PN L (ω) = χ(1) + 2χ(2) E(0) E(ω).

(2.33)

La ecuaci´on anterior significa que el indice de refracci´on de afecta por la presencia del campo el´ectrico constante n2 = 1 + 4πχef f = 1 + 4πχ(1) + 4π(2χ(2) E(0)) = (n + ∆n)2 .

(2.34)

Usando el hecho de que los efectos no lineales son peque˜ nos (∆n ¿ n) tenemos que ∆n2 ' 2n∆n = 8πχ(2) E(0),

(2.35)

luego ∆n =

2.2

4πχ(2) E(0) n

(2.36)

Mezcla de dos ondas

Consideremos dos ondas monocrom´aticas incidentes E = E(ω1 )eiω1 t + E(ω2 )eiω2 t + c.c.

(2.37)

£ PN L = χ(2) E(ω1 )2 ei2ω1 t + E(ω2 )2 ei2ω2 t + 2E(ω1 )E(ω2 )ei(ω1 +ω2 )t ¤ +2E(ω1 )E(ω2 )∗ ei(ω1 −ω2 )t + c.c.

(2.38)

Entonces

De aqui se deducen las polarizaciones no lineales PN L (2ω1 ) PN L (2ω2 ) PN L (ω1 + ω2 ) PN L (ω1 − ω2 )

= = = =

χ(2) E(ω1 )2 χ(2) E(ω2 )2 2χ(2) E(ω1 )E(ω2 ) 2χ(2) E(ω1 )E(ω2 )∗

(2.39) (2.40) (2.41) (2.42)

´ CHAPTER 2. ELEMENTOS DE OPTICA NO LINEAL

66

Para que estos efectos sean observables es requisito que se cumpla la condici´on de phase matching . Consideremos la mezcla de dos ondas planas ~

(2.43)

~

(2.44)

E(ω1 ) = A1 e−ik1 ·~r E(ω2 ) = A2 e−ik2 ·~r Entonces

~

~

PN L (ω1 + ω2 ) = 2χ(2) A1 A2 e−i(k1 +k2 )·~r .

(2.45)

Las leyes de conservaci´on de la energia y del momentum imponen que ω3 = ω1 + ω2 ~k3 = ~k1 + ~k2 .

(2.46) (2.47)

La ley de dispersi´on da una relaci´on entre los vectores de onda y las frecuencias ki =

2.2.1

ni ωi . c

(2.48)

Phase matching en problemas unidimensionales

Las leyes de conservaci´on toman la forma ω3 = ω1 + ω2 n3 ω3 = n1 ω1 + n2 ω2

(2.49) (2.50)

De las ecuaciones anteriores se obtiene (n3 − n1 )ω1 = (n2 − n3 )ω2 .

(2.51)

En la mayor parte del espectro visible en medios transparentes se cumple que dn/dω > 0. Esto implica que si ω1 < ω2 < ω3 entonces n1 < n2 < n3 y la ecuaci´on (2.51) no tiene soluci´on. La soluci´on es posible, y ello determina la observabilidad del efecto de suma de dos ondas, en regiones espectrales de dispersi´on an´omala (dn/dω < 0), las que existen cerca de las bandas de absorci´on.

Amplificador

Up converter

ω1

ω3

ω2

ω1 ω3

ω1

Figure 2.3: Dispositivos que utilizan phase matching.

2.2.2

Aplicaciones

En el seno del material optico estan presentes las tres frecuencias ω1 , ω2 y ω3 , de modo que son posibles todas las combinaciones entre estas. En la Figura 2.3 se ilustran dos aplicaciones.

´ 2.3. FENOMENOS DE TERCER ORDEN

2.3

67

Fen´ omenos de tercer orden

En medios con simetria de inversi´on se cumple que χ(2) = 0, por lo cual se hace necesario considerar los procesos de orden 3, los cuales est´an regidos por la relaci´on (3)

PN L = χ(3) E 3 .

2.3.1

(2.52)

Efecto electro´ optico

Si tenemos un campo constante intenso y una onda monocrom´atica, dado por (2.37), con E(0) À E(ω) obtenemos ¡ ¢ (3) PN L (ω) = χ(3) 3E(0)2 + 2|E(ω)|2 E(ω) ' 3E(0)2 χ(3) E(ω). (2.53) El resultado anterior significa una modificaci´on al indice de refracci´on lineal, dada por ∆n =

2.3.2

6πχ(3) E(0)2 n

(2.54)

Tercer arm´ onico y efecto Kerr ´ optico

Si no hay campo electrico constante, solamente la onda incidente de frecuencia ω, entonces se obtiene (3)

PN L (3ω) = χ(3) E(ω)3

Generaci´on del tercer arm´onico.

(3)

PN L (ω) = 3χ(3) |E|2 E(ω) Efecto Kerr.

(2.55) (2.56)

El resultado anterior implica un cambio del ´ındice de refracci´on ∆n =

6πχ(3) |E(ω)|2 = n2 I n

(I ∝ |E(ω)|2 ).

(2.57)

∆ϕ

χ(3) Figure 2.4: Esquema de dispositivo controlado por el efecto Kerr. Veamos como influye el efecto Kerr en una onda plana. La fase acumulada durante un recorrido L es 2π (n + n2 I)L. (2.58) ϕ(L) = λ

´ CHAPTER 2. ELEMENTOS DE OPTICA NO LINEAL

68

Esto implica que el cambio de fase asociado a la intensidad de la luz es ∆ϕ(L) =

2π n2 IL. λ

(2.59)

Esto permite realizar dispositivos de interferencia que regulen la intensidad de la luz emitida mediante un cambio de fase, que es a su vez dictado por la intensidad de la luz incidente. Esto se ilustra en la figura 2.4. En la Figura 2.4, el dispositivo sealizado por χ(3) , puede ser una fibra optica con los siguientes parametros L = 1 m, A = 10−2 mm2 , n2 = 10−10 cm2 /W, λ = 1µmm.

(2.60)

Segun la ecuaci´on (2.59), la potencia reeuqrida P = IA para provocar una interferencia destructiva ∆ϕ = π es P = 1.5W .

2.3.3

Autoenfoque

n=n0 + n2 I

Otra aplicaci´on importante es el autoenfoque, ilustrada en la Figura 2.5. Un haz de luz se describe generalmente por un perfil de intensidad gausiano. Al atravesar un medio cuyo indice de refracci´on esta dado por N = n0 + n2 I, los frentes de onda se encorvan, de modo que el producto n2 IL sea constante en el frente de onda (vease (2.59)) y al salir del medio los rayos convergen en un punto focal.

I

Figure 2.5: Efecto de autoenfoque de un haz de luz de perfil gausiano. Demostremos esto anal´ıticamente. Consid´erese la Figura 2.6 y el angulo θ(y) formado por un rayo con el eje Z. Considerendo la refracci´on en una lamina paralela a Z de ancho dy, la ley de Snell se lee n(y) cos θ(y) = n(y + dy) cos θ(y + dy).

(2.61)

esto lleva a la ecuaci´on diferencial dθ dn = n tan θ . dy dy

(2.62)

´ 2.3. FENOMENOS DE TERCER ORDEN

69

y

y+dy y

θ (y+dy) n(y) θ (y)

z

Figure 2.6: Refracci´on de un haz de luz en un medio de indice de refracc´on variable.

Ahora bien, tan θ = dy/dz, siendo y(z) la trayectoria del rayo de luz. Usaremos la aproximaci´ on paraxial, o sea, que consideramos rayos casi paralelos al eje focal (θ ∼ 0). En este caso tan θ ' θ y la ecuaci´on diferencial toma la forma 1 dn d2 y = 2 n dy dz

(2.63)

El haz de luz tiene el perfil I = I0 exp(−αy 2 ). Si consideramos la regi´on cercana al eje focal, podemos aproximarla por I = I0 (1 − αy 2 ). El indice de refracci´on queda igual a µ ¶ β 2y2 n2 n = n0 + n2 I = n0 1 − , β = αI0 . (2.64) 2 n0 Introduciendo la formula anterior en la ecuaci´on en (2.63) se obtiene y 00 + β 2 y = 0

(2.65)

Y por tanto la soluci´on es

θ0 sin(βz) (2.66) β En una fibra ´ optica la longitud L es mucho mayor que y(z), de modo que los rayos de luz siguen una trayectoria sinusoidal. y0 es la altura a la que penetra el rayo, y θ0 es el angulo que forma al penetrar. En una lente delgada ocurre lo contrario, el grosor de la lente es L ¿ y(z). y(z) = y0 cos(βz) +

y(L) ' y0 θ(L) = βy0 sin βL ' y0 β 2 L tan θ(L) = dz

(2.67) (2.68)

Despues de salir de la lente, el rayo de luz viaja en linea recta. Sea f el punto en que y(f ) = 0 f=

y0 1 n0 = 2 = . tan θ(L) β d n2 αLI0

(2.69)

f es independiente de la altura y0 con que inciden los rayos paralelos, lo que implica que f es una distancia focal.

70

´ CHAPTER 2. ELEMENTOS DE OPTICA NO LINEAL

Chapter 3 Elementos de biof´ısica 3.1

Biomec´ anica

En esta secci´on aplicaremos las leyes de la mec´anica para estimar las fuerzas que actu´ uan sobre algunos musculos y huesos del cuerpo humando.

3.1.1

Fuerzas que act´ uan sobre el f´ emur y la cadera

La cadera esta formado por un conjunto de huesos que en la infancia est´an unidos por cart´ılagos y en la adultez quedan soldados. La Fig. 3.1 ilustra este caso. La cadera se puede dividir en tres partes: el ilio, el isquion (este es el hueso que se apoya en posici´on sentada) y el pubis. Las partes derecha e izquierda se unen por detras con el sacro y el coxis, que son las partes inferiores de la columna vertebral, mientras que por delante se unen en el pubis. El conjunto completo se denomina pelvis.

Iliac crest

Pubic crest

Ischial spine

Figure 3.1: Vista frontal y lateral de los huesos de la cadera. Prestaremos especial atenci´on al acetabulo, que es la cavidad en que se ajusta la cabeza del f´emur (hueso del trutro :-) ) y cuya forma redondeada permite el movimiento articulado de la pierna (Vease la Fig. 3.2). 71

72

CHAPTER 3. ELEMENTOS DE BIOF´ISICA

Figure 3.2: Vistas frontal (izq) y trasera (cen) del f´emur. Estructura de la cabeza del f´emur (der). Cubriendo la cabeza del femur se encuentra la epifisis, que es solo la parte superficial del hueso. En casos patologicos (epifisiolisis) la epifisis se puede desprender parcialmente del femur si est´a sujeta a esfuerzos de corte (o cizalla). El f´emur tiene varias protuberancias que sirven de puntos de inserci´on a los m´ usculos, de las cuales en la Fig. 3.2 se destaca el trocanter mayor. En el trocanter mayor se insertan los tendones que de cinco m´ usculos que 1 colaboran en el movimiento de abducci´on. Los m´ usculos abductores m´as importantes son el gl´ uteo medio y el gl´ uteo minimo. Los extremos opuestos se pegan en varias partes del ilio. Hay otros musculos presentes en esta zona, cuya funci´on es rotar la pelvis y no desempe˜ nan ning´ un rol en la abducci´on. El tercer m´ usculo importante es el llamado tensor de la fascia lata, que se une al femur ligeramente bajo el trocanter mayor y por el otro extremo se une a la fascia lata y el tracto iliotibial, cerca de la rodilla. La importancia relativa de estos m´ usculos y su linea de acci´on efectiva ha sido estudiada 2 mediante radiograf´ıas. Se ha establecido que la l´ınea de acci´on de de los m´ usculos abductores ◦ pasa por el trocanter mayor y forma un ´angulo de 70 respecto a la horizontal (vea Fig. 3.3).

3.1.2

Fuerzas en equlibrio sobre un pie

Con los datos anteriores se pueden calcular las fuerzas realizadas sobre el trocanter mayor y sobre la cabeza del femur. Consideremos el caso de una persona en equilibrio sobre un solo pie. El esquema de fuerzas se ilustra en la Fig. 3.3. N es la fuerza de reacci´on del suelo y es igual al peso del cuerpo W . F1 es la fuerza neta de los musculos abductores, R es la fuerza de reacci´on que ejerce el acet´abulo sobre la cabeza del f´emur, y WL es el peso de la 1 La abducci´ on es el movimiento de abrir los muslos lateralmente, y algunos creen que los extraterrestres le hacen eso a los humanos que llevan a sus naves. La aducci´ on es el movimiento contrario, o sea, cerrar los muslos. 2 Vease G. B. Benedek, Physics with illustrative examples from medicine and biology (New York : SpringerVerlag, 2000), Capitulo 3, y referencias citadas ah´ı.

´ 3.1. BIOMECANICA

73

Figure 3.3: Esquema de las fuerzas que actuan sobre la pierna cuando el apoyo es sobre un solo pie. propia pierna, actuando sobre el centro de masa, y que se estima WL ' 17 W . Apliquemos las condiciones de equilibrio est´atico a la pierna 1 F1 sin 70◦ − Ry − W + W = 0, (3.1) 7 F1 cos 70◦ − Rx = 0, (3.2) 1 (F1 sin 70◦ )(2.75) + ( W )(1.25) − W (4.25) = 0 (torque respecto al acetabulo). (3.3) 7 Resolviendo el sistema de ecuaciones se obtiene F1 = 1.58W, Rx = 0.54W, Ry = 2.34W, R =

q

(3.4) Rx2 + Ry2 = 2.4W, φ = 13◦ .

(3.5)

As´ı encontramos que sobre la cabeza del f´emur, en el instante en que se est´a en equilibrio sobre un solo pie, act´ ua una fuerza igual a 2.4 veces el peso del cuerpo. El resultado anterior tiene implicaciones clinicas y anat´omicas. La Fig. 3.4 muestra una imagen de rayos X y un esquema de la microestructura del cuello femoral. La estructura osea forma una red esponjosa (trab´ecula) cuya direcci´on de crecimiento coincide con la l´ınea de acci´on de la fuerza R. Cuando la pelvis rota hacia arriba y hacia abajo en el plano vertical, resulta que R se mantiene dirigida a lo largo de la trab´ecula. De este modo la evoluci´on se ha dirigido a evitar grandes esfuerzos de cizalla sobre el cuello del f´emur y sobre la epifisis.

CHAPTER 3. ELEMENTOS DE BIOF´ISICA

74

Figure 3.4: Estructura osea del femur e imagen de rayos X. En el caso de haber lesiones en los m´ usculos abductores ocurren cambios interesantes. Es claro que si se reduce la fuerza F1 que pueden ejercer los musculos es necesario variar la posici´on del cuerpo para mantener el equilbrio. El enfermo corrige esto instintivamente (para evitar dolor y para mantener el equilibrio) moviendo su centro de gravedad hacia el lado del musculo afectado de modo que se reduzca el torque de la normal N . Esto es lo que se llama marcha antalgica, o sea, una marcha encorvada para evitar el dolor. Si el encorvamiento es tal que la linea de accion del peso y la normal pasan por la cabeza del femur y a lo largo de la pierna, entonces el brazo de las fuerzas N y WL se hace 0 y las ecuaciones de equilibrio se satisfacen con F1 =0. El precio que se paga es que apareen esfuerzos de corte en el cuello femoral y el consecuente riesgo de fractura. Si la afecci´on en los m´ usculos abductores ocurre en etapas tempranas de la vida y se prolonga, entonces el crecimiento del cuello del femur se orienta hacia arriba en respuesta a la fuerza vertical R. Esto se denomina coxa valga. Como resultado, un femur se hace mas largo que el otro, la pelvis queda en posicion inclinada, y se producen deformaciones en la columna vertebral.

3.1.3

Efecto de un bast´ on

Como hemos visto, la fuerza ejercida por los musculos abductores para mantener el equilibrio depende notablemente del torque ejercido por la fuerza de reacci´on del suelo sobre el pie. Si se reduce el brazo de esta fuerza, se reduce la fuerza F1 . Consideremos un hombre que tiene afectado los musculos abductores del lado derecho. Durante la caminata, al levantar el pie izquierdo se da el caso calculado en la seccion previa, y si hay afectaciones en los musculos abductores derechos, necesitar´a inclinarse hacia la derecha para mantener su momentaneo equilibrio. Una alternativa es usar un bast´on en la mano izquierda (n´otese que es el lado opuesto al musculo afectado), de modo que reduzca la magnitud y el brazo de la fuerza N . La Fig. 3.5 muestra el esquema de fuerzas para el hombre como un todo y para la pierna

´ 3.1. BIOMECANICA

75

Figure 3.5: Esquema de las fuerzas que actuan sobre la pierna cuando el apoyo es un solo pie y un baston. afectada. Llamamos C a la fuerza trasmitida al cuerpo por el baston. Un hombre que pese 200 libras puede ejercer cerca de 30 lb sobre el baston sin resultar un esfuerzo excesivo, lo cual representa C ' 16 W . Las otras fuerzas externas son el peso W y la reacci´on del suelo N . Ademas, podemos estimar en 12 pulgadas la distancia horizontal del baston a la columna vertebral. (Fig. 3.5 central). En estas condiciones podemos calcular la posici´on del pie y la fuerza N necesaria para mantener el equilibrio. Llamemos l a la distancia horizontal del punto de apoyo respecto a la linea media del cuerpo. Las condiciones de equilibrio para el cuerpo completo son W − W = 0 (fuerza neta), 6 W N l − (1200 ) = 0 (torque respecto al centro de masa). 6 N+

(3.6) (3.7)

Las soluciones son

5 N = W, l = 2.400 . (3.8) 6 Con estos datos podemos resolver las ecuaciones de equilibrio para la pierna (Fig. 3.5 derecha), con lo cual se obtiene F1 = 0.60W,

R = 1.26W,

φ = 9◦ .

(3.9)

El bast´on ha permitido reducir la fuerza sobre la cabeza femoral de R = 2.34W a R = 1.26W (En un hombre de 200 lb esto es pasar de 468 lb a 252 lb) gracias al m´odico esfuerzo de 30

76

CHAPTER 3. ELEMENTOS DE BIOF´ISICA

lb en el brazo opuesto. Este beneficio desempe˜ na un rol importante en la rehabilitaci´on de pacientes que han sufrido cirug´ıa de la cadera o fractura de femur.

3.2

Fuerzas que act´ uan sobre las v´ ertebras lumbares

La Fig. 3.6 muestra un diagrama de la columna vertebral, donde el frente est´a hacia la derecha y la espalda hacia la izquierda. La parte inferior es el coxis, seguida del sacro. Encima se encuentran cinco vertebras lumbares, numeradas desde arriba hacia abajo (la m´as cercana al sacro es la quinta, la mas lejana es la primera). Las cinco v´ertebras forman una curva llamada lordosis lumbar. Encima se hallan doce vertebras tor´axicas, dispuestas con una curvatura inversa denominada cifosis. Encima se encuentran las siete vertebras cervicales. El sacro est´a soldado rigidamente a la pelvis en las personas adultas, y su geometria se describe por el angulo lumbosacral (ver Fig). Este angulo determina la magnitud de la lordosis lumbar. La pelvis puede girarse hacia adelante si los musculos abdominales estan debilitados, lo cual incrementa la lordosis. El angulo lumbosacral normal es cercano a 30◦ .

Figure 3.6: Izquierda: esquema de la columna vertebral humana. Derecha: vista lateral de la de la cadera y la parte inferior de la columna vertebral. Examinemos una v´ertebra. La Fig.3.7 en su parte izquierdaa muestra dos vertebras vistas desde el lado. La parte anterior es la que soporta el peso del cuerpo. Para esto se ayuda de unos discos que est´an entre cada dos vertebras consecutivas. Este disco es un sistema fluido autocontenido (nucleo pulposo) que absorbe choques y trasmite la presion uniformemente, permite deformacion del espacio intervertebral y por tanto permite el movimiento. Las partes superior e inferior del disco intertertebral son cartilagos y estan unidos al hueso de las vertebras. La pared lateral, llamada anillo fibroso, se compone de una docena de capas de fibras de colageno. El interior del disco es un gel viscoelastico, 80% agua. Durante la segunda

´ ´ 3.2. FUERZAS QUE ACTUAN SOBRE LAS VERTEBRAS LUMBARES

77

Figure 3.7: Esquema de una v´ertebra. decada de vida los discos intervertebrales tienen suministro de sangre, y luego se alimentan por difusion de linfa. La elasticidad es proporcionada en buena medida por el anillo fibroso, la cual degenera con la edad o por lesiones producidas por sobrecargas. La hernia discal consiste en la hernia del anillo fibroso y la extrusion del gel, presionando sobre los nervios vecinos o la medula espinal y produciendo dolor. El estrechamiento del espacio intervertebral termina provocando da˜ no en las vertebras por contacto entre ellas. La Fig. 3.8 muestra la curva de elasticidad de los discos intervertebrales lumbares. Con una carga mayor de 1500 kg, el disco se rompe. Para las vertebras superiores la carga maxima disminuye, siendo de 320 kg para las vertebras cervicales. El area de los discos intervertebrales decrece, de modo que la presion maxima es casi constante, cercana a 1.1 kg/mm2 . En la parte derecha de la

Figure 3.8: Carga soportada por una disco lumbar vertebral en funci´on de la deformaci´on, para personas de 40 a 59 a˜ nos.

CHAPTER 3. ELEMENTOS DE BIOF´ISICA

78

Figure 3.9: Lo que ocurre cuando hay degeneracion del disco intervertebral. Fig.3.7 se ve una vertebra desde arriba. Las protuberancias (denominadas procesos), actuan como puntos de inserci´on para los musculos. El desplazamiento horizontal de las vertebras es limitado por la estructura osea, los ligamentos y los musculos. Observese el ligamento longitudinal posterior, en la vista superior. Este ligamento impide impide la extrusion del fluido del disco intervertebral hacia la cavidad, que es por donde pasa la medula espinal y por donde salen los nervios. Resulta que este ligamento es mas delgado en las vertebras lumbares, en especial en la quinta, precisamente donde mas necesario es. Esto es obviamente un punto debil de la anatomia humana, un rezago de la cuadrupedia que la evoluci´on no ha alcanzado a corregir en el mill´on de a˜ nos de posici´on erecta. La Figura 3.9b muestra lo que ocurre cuando hay degeneraci´on lumbar. El disco intervertebral presiona la raiz del nervio y produce dolor. Es por eso que las lesiones ocasionadas por esfuerzos al cargar pesos se producen generalmente en la region lumbar. Examinemos las fuerzas que actuan entre las v´ertebras cuando el cuerpo se alza un peso. En la Fig. 3.10 se muestra el diagrama de fuerzas que act´ uan sobre la columna al levantar el tronco. W1 = 0.4W es el peso del torax, W2 es el peso de los brazos, cabeza y carga sujeta en las manos. En ausencia de carga puede W2 = 0.2W . R es la fuerza de reacci´on ejercida por el sacro sobre el disco lumbo-sacral. La fuerza Fe es la resultante de la tensi´on ejercida por los musculos sacroespinal y erector espinal, que se insertan en el ilio y bajo sacro por un extremo y en cuatro de las vertebras toraxicas por arriba. Estos son, obviamente, los que hacen el trabajo de elevaci´on. La geometria de todos esos musculos ha sido estudiada mediante radiografias y puede ser reemplazado por la fuerza Fe , en el punto de inserci´on indicado y formando un angulo de 12◦ con la columna. Considerando un angulo θ = 30◦ y aplicando las ecuaciones del equilibrio estatico se obtiene Fe = 2.5W,

R = 2.74W,

φR = 35◦ .

(3.10)

Como se ve, son fuerzas considerables. Para un hombre de 200 lb, la fuerza de compresi´on es R = 540 lb, equivalente a una carga de 250 kg. La Fig. 3.8 informa que esto corresponde a una compresi´on del 20 % del espacio intervertebral. Considerando constante el volumen, esto implica un aumento del 10 % en el radio del disco, que va a presionar en la medula espinal.

´ ´ 3.2. FUERZAS QUE ACTUAN SOBRE LAS VERTEBRAS LUMBARES

79

Figure 3.10: Esquema de la acci´on de los m´ usculos de la espalda. Si con las manos se toma un peso de un quinto del peso corporal (40 lb), entonces W2 = 0.2W + 0.2W , y se obtiene una fuerta de compresi´on igual a R = 4.07W (equivalente a una carga de 370 kg en el ejemplo visto). La fuerza efectuada por los musculos erectores es Fe = 3.74W . Las consecuencias de esto son claras. Con esto se puede entender facilmente porque las posiciones correcta e incorrecta de cargar un peso son las mostradas en la Fig. 3.11.

Figure 3.11: Derecha: posicion correcta de cargar un peso. Izquierda: posicion incorrecta.

80

CHAPTER 3. ELEMENTOS DE BIOF´ISICA

Chapter 4 F´ısica de part´ıculas elementales

81

82

CHAPTER 4. F´ISICA DE PART´ICULAS ELEMENTALES

Chapter 5 F´ısica del n´ ucleo at´ omico

83

84

´ ´ CHAPTER 5. F´ISICA DEL NUCLEO ATOMICO

Chapter 6 Mapas discretos En esta parte estudiaremos la evoluci´on de las relaciones de recurrencia del tipo xn+1 = M (xn )

(6.1)

A este tipo de relaci´on se le llama mapa. Un ejemplo cl´asico, que utilizaremos para ilustrar los conceptos, es el mapa log´ıstico, definido por M (x) = rx(1 − x),

r > 0.

(6.2)

El mapa log´ıstico permite modelar simplificadamente la evoluci´on de la poblaci´on de una especie cuyo suministro alimentario presenta un tope. Como veremos, este mapa presenta caos y orden, en dependencia del valor del par´aametro r. El mapa logistico es unidimensional, pero en general xn puede representar un punto en un espacio de mas dimensiones. Muchos problemas de la f´ısica son gobernados por ecuaciones diferenciales ordinarias, cuya soluci´on num´erica implica discretizar el tiempo segun tn = n∆t. Si como resultado se obtienen relaciones del tipo x(tn+1 ) = f [x(tn )], entonces estamos nuevamente en presencia de un mapa.

6.1

El mapa log´ıstico

Estudiemos gr´aficamente la evoluci´on del mapa log´ıstico. La Fig. 6.1 muestra la secuencia de valores que toma xn a partir de una condicion inicial dada. En la parte izquierda la condicion inicial es x0 = 0.5. La linea vertical que parte de (0.5,0) hasta (0.5,M(0.5)) [es decir, (x0 , x1 )] indica la evaluaci´on del mapa por primera vez. Luego se traza una linea horizantal que intersecta la recta y = x, en el punto (x1 , x1 ). Trazando la vertical hasta la funcion M (x) se obriene el punto (x1 , M (x1 ) = (x1 , x2 )]. Continuando el proceso obtenemos la secuencia de puntos xn , xn+1 . De la figura se deduce que la sucesion de valores xn converge hacia el valor 0 si x0 = 0.1 y r = 0.9, y converge hacia un cierto valor de x > 0.5 en el caso de que r = 2.5. En ambos casos, se obtiene siempre el mismo punto de convergencia si se cumple que la condici´on inicial es 0 < x0 < 1. El conjunto de puntos {x1 , x2 , x3 , ....} se denomna trayectoria. El conjunto al que convergen asintoticamentes las trayectorias es llamado atractor. En el caso r = 0.9 el atractor es el punto x = 0, para 85

86

CHAPTER 6. MAPAS DISCRETOS

M(x)

M(x) 1

r=0.9

0.25

r=2.5

0.8

0.2 0.6 0.15 0.4 0.1 0.2

0.05 0.2

0.4

0.6

1

0.8

x

0.2

0.4

0.6

0.8

1

x

Figure 6.1: Evoluci´on del mapa log´ıstico a partir de una condici´on inicial. Derecha: r = 0.9, x0 = 0.5. Izquierda: r = 2.5, x0 = 0.1.

x 1 0.8 0.6 0.4 0.2

1

3

2

4 x 1

r

0.8 0.6 0.4 0.2

3.5

3.6

3.7

3.8

3.9

r

Figure 6.2: Diagrama de bifurcaciones para el mapa log´ıstico, si se cumple que la condicion inicial es x0 ∈ (0, 1).

6.1. EL MAPA LOG´ISTICO

87

r = 2.5 el atractor es el punto x que se ve en la figura. Se puede tener una idea global de los atractores en funci´on del valor del par´ametro si se hace un diagrama de bifurcaciones. La Fig. 6.2 muestra el atractor para cada valor del parametro r, siempre que 0 < x0 < 1. El atractor es x = 0 si 0 < r < 1, es un punto 0 < x < 0.8 si 1 < r < 3, dos puntos si 3 < r < 3.45, luego 4, etc, hasta llegar a formar una banda continua en cierto rango de r. En esta banda estamos en precencia de caos. Como indica la trayectoria, se puede ir a cualquier parte. El grafico es realmente rico, en la parte de abajo se ve una ampliaci´on en el rango 3.5 < r < 4, en la que se nota que las bifurcaciones van en aumento al acercarse al regimen de caos, y la aparicion de ventanas de orden dentro del caos. La Fig. 6.2 se genera aplicando el siguiente algoritmo para cada valor del parametro r • Condici´on inicial x0 = 0.1. • Se aplico el mapa 200 veces, o sea, se genera la trayectoria hasta x200 . • Se aplico el mapa 200 veces mas y se graficaron los puntos (r, xn ). El segundo paso del algoritmo consiste en eliminar el transiente, de modo que lo puntos graficados, generados en el tercer paso, ya corresponden a la tendencia asint´otica de las trayectorias, o sea, el atractor. La figura generada no depende del valor preciso de x0 , siempre que x0 ∈ (0, 1). El intervalo (0, 1) es en este caso la cuenca1 de atracci´ on. Puede comprobarse siguiendo las trayectorias gr´aficamente que si x0 < 0 o x0 > 1, entonces el atractor es −∞. Veamos como puede entenderse el resultado anterior usando procedimienos anal´ıticos. Si la trayectoria converge hasta cierto punto x∗ = lim xn

(6.3)

x∗ = M (x∗ ).

(6.4)

n

se cumple que Los puntos que satisfacen la ecuaci´on anterior se denominan puntos fijos del mapa M . Para el mapa log´ıstico las soluciones de la ecuaci´on 6.4 son x(1) = 0,

1 x(2) = 1 − . r

(6.5)

Los puntos fijos no son necesariamente atractores. Examinemos como evoluciona una trayectoria cerca de un punto fijo. Sea xn = x∗ + δn ¯ ¯ dM ¯ x∗ + δn+1 = xn+1 ' M (x∗ + δn ) = M (x∗ ) + δn . (6.6) dx ¯x=x∗ Utilizando 6.4 obtenemos la relaci´on ¯ ¯ ¯ ¯ ¯ δn+1 ¯ ¯ dM ∗ ¯ ¯ ¯ ¯ ¯ ¯ δn ¯ = ¯ dx (x )¯ 1

En ingl´es, basin, que significa lavamanos.

(6.7)

88

CHAPTER 6. MAPAS DISCRETOS M(x) 1

M2 (x) 1

r=3.2

0.8

0.8

0.6

0.6

0.4

0.4

0.2

0.2

0.2

0.4

0.6

1

0.8

x

r=3.2

0.2

0.4

0.6

0.8

1

x

Figure 6.3: Orbita estable de periodo 2 n el mapa logistico. De aqui se deduce que δn (y por tanto xn −→ x∗ ) tiende a cero solamente si |M 0 (x∗ )| < 1.

(6.8)

En este caso se dice que el punto fijo es estable y es un atractor del mapa. Evaluando esta condici´on en los puntos fijos obtenidos arriba, obtenemos M 0 (x(1) ) = r,

M 0 (x(2) ) = 2 − r.

(6.9)

De aqui se deduce que el punto fijo x(1) es estable en el rango r ∈ (0, 1) y x(2) es estable en r ∈ (1, 3). Estos son los valores que se observan en el diagrama de bifurcaci´on de la Fig. 6.2. Otro asunto ocurre para r > 3. La Fig. 6.3 ilustra una trayectoria para r = 3.2. Se aprecia que despu´es de un transiente la trayectoria se estabiliza en los valores alternantes 0.513 y 0.799. Este tipo de atractor de llama ´ orbita de per´ıodo 2, y se generaliza a periodo N de manera directa. Los dos valores de esta orbita son los que se aprecian en el diagrama de bifurcaciones. A partir de cierto valor de r ∼ 3.44 las orbitas son de periodo 4 y al aumentar r se siguen multiplicando hasta llegar al caos. Entendamos una ´orbita de per´ıodo 2. La alternancia de la trayectoria entre dos puntos x(1) y x(2) significa que ellos son puntos fijos estables del mapa compuesto M 2 (x) = M (M (x)). Es decir M 2 (x∗ ) = x∗ . (6.10) La Fig. 6.3 muestra una trayectoria del mapa M 2 que converge a uno de los puntos fijos estables. N´otese que la ecuaci´on 6.10 incluye las soluciones de 6.4. Esto se comprueba facilmente sustituyendo 6.4 en 6.10. Comparando las partes izquierda y derecha de la Fig. 6.3 que efectivamente las soluciones de la parte izquierda (dadas por la interseccion del mapa con la recta y = x) estan presentes en la figura derecha, aunque son inestables. La condici´on de estabilidad para una ´orbita de periodo 2 que alterna entre los puntos x(1) y x(2) esta dada por |dM 2 /dx| < 1. Aplicando la regla de la cadena se obtiene ¯ ¯ ¯ ¯ ¯ ¯ dM dM d ¯ ¯ . = (6.11) M (M (x))¯¯ dx dx ¯ (2) dx ¯ (1) (1) (1) x=x

x

=M (x

)

La generalizaci´on para orbitas de periodo m es obvia ¯ ¯ m Y dM ¯¯ d m ¯¯ = M (x)¯ ¯ (k) . dx dx x=x(1) x=x k=1

x

(6.12)

6.1. EL MAPA LOG´ISTICO

89

λ 1 0.5

1

1.5

2

2.5

3

3.5

4

r

-1 -2 -3 -4 λ 1 0.5 3.4

3.5

3.6

3.7

3.8

3.9 -0.5

r

-1 -1.5 -2

Figure 6.4: Exponente de Lyapunov del mapa log´ıstico, calculado mediante la ecuaci´on 6.16, aproximando m = 200. Esto implica que la desviaci´on δn = xn − x(1) cumple ¯m ¯ ¯Y dM ¯¯ ¯ ¯ ¯ ¯ |δn+m | = ¯ ¯ |δn | . ¯ ¯ dx (k) ¯ x=x k=1 Los valores |δn | admiten la soluci´on

|δn | = eλn ,

que sustituida en la condici´on 6.13 implica que ¯ ¯ m 1 X ¯¯ dM (k) ¯¯ (x )¯ . ln λ= m k=1 ¯ dx

(6.13)

(6.14)

(6.15)

En la pr´actica, los x(k) se pueden remplazar por los puntos consecutivos de la trayectoria una vez eliminado el transiente. El numero λ se llama exponente de Lyapunov. Si λ > 0 hay caos, en caso contrario hay orbitas m-periodicas estables. El caso de las trayectorias ca´oticas se puede considerar como orbitas de periodo infinito, con un exponente de Lyapunov igual a ¯ ¯ m 1 X ¯¯ dM (k) ¯¯ (6.16) (x )¯ . λ = lim ln ¯ m→∞ m dx k=1

90

CHAPTER 6. MAPAS DISCRETOS

La ecuaci´on anterior (6.16) incluye a 6.15. El exponente de Lyapunov del mapa log´ıstico se muestra en la Fig. 6.4. Cuando λ < 0 estamos en presencia de puntos fijos o orbitas de per´ıodo n. Cuando λ = 0 tenemos una bifurcaci´on y generalmente se acompa˜ na de un cambio de pendiente producto de que cambia el atractor. λ → −∞ corresponde puntos superestables, o sea, que M 0 (x∗ ) = 0, y λ > 0 indica caos. Este grafico revela limpiamente la existencia de caos o orden. En la ampliacion para r ∈ (3.4, 4.0) (comp´arese con Fig. 6.2)se observa la existencia de muchas ventanas de orden dentro del caos.2 El exponente de Lyapunov nos da otra informaci´on. Si tenemos dos trayectorias que parten desde cerca del atractor, tales que |x0 − y0 | = δ ¿ 1, entonces |xn − yn | = δenλ .

(6.17)

Es decir, si λ > 0 las trayectorias divergen de forma exponencial hasta que se alcanza cierto l´ımite determinado por la imagen del mapa.3

2 3

Un estudio detallado revela que es infinito el numero de ventanas de orden dentro del caos. En el mapa logistico |xn − yn | < 1.

Chapter 7 Caos en ecuaciones diferenciales ordinarias En esta parte se pretende extender los conceptos de mapas discretos a los sistemas de ecuaciones diferenciales ordinarias. Consideremos un sistema de ecuaciones diferenciales de primer orden, escrito genericamente en la forma dx = F (x, p), x(0) = x0 , (7.1) dt donde p es un conjunto de parametros y x es un vector de N dimensiones, cuyo conjunto de valores posibles es llamado espacio de fase. Puede demostrarse que si las funciones F y dF/dt son continuas en una vecindad S de x0 , entonces existe algun intervalo (−τ, tau) de medida no nula en el cual el sistema de ecuaciones tiene soluci´on unica. Si la soluci´on es u ´nica, entonces dos trayectorias x(t) no se cruzan nunca, siempre que la funci´on F (x, p) no dependa explicitamene de t. El sitema de ecuaciones se puede resolver num´ericamente usando el algoritmo de Runge-Kutta. dx dt k1 k2 k3 k4

= f (x, t)

(7.2)

= = = =

(7.3) (7.4) (7.5) (7.6)

∆tf (x, t) ∆tf (x + k1 , t + ∆t/2) ∆tf (x + k2 /2, t + ∆t/2) ∆tf (x + k3 , t + ∆t) k1 + 2k2 + 2k3 + k4 x(t + ∆t) = x(t) + + O((∆t)6 ) 6

(7.7)

Un ejemplo cl´asico es el sistema de Lorenz, obtenido en 1963 en un estudio de fenomenos de convecci´on. dx = σ(y − x) dt dy = −xz + rx − y dt dz = xy − bz. dt 91

(7.8) (7.9) (7.10)

92

CHAPTER 7. CAOS EN ECUACIONES DIFERENCIALES ORDINARIAS

y

20 10

0 -10 -20

40 30 z 20 10 0 -10 0 x

10

Figure 7.1: Trayectoria caotica del sistema de Lorenz. σ = 10, b = 8/3, r = 28, t ∈ (60, 100). Este sistema presenta caos para ciertas condiciones iniciales y valores de los par´ametros, por ejemplo, σ = 10, b = 8/3, r = 28. En la Fig. 7.1 se muestra parte del attractor ca´otico del sistema de Lorenz. El caracter ca´otico se pone de manifiesto considerando dos trayectorias ~x1 (t) y ~x2 (t) tales que δ(t0 ) = |~x1 (t0 ) − ~x2 (t0 )| = 10−10

(7.11)

con t0 suficientemente largo para eliminar el transiente y estar ya en el atractor. Si se plotea lg δ(t) vs t, con t0 = 55 se obtiene el grafico de la Fig. 7.2. lg δ -2

-4

-6

60

65

70

75

t

-10

Figure 7.2: Evolucion de la diferencia de dos trayectorias con condiciones iniciales casi iguales en el sistema de Lorenz. La linea recta es un ajuste a la relacion lg δ(t) = λt + C.

93 Definiendo el exponente de Lyapunov mediante la relaci´on δ(t) = δ(t0 )eλ(t−t0 ) ,

(7.12)

se obtiene para el caso de la Fig. 7.2 el valor λ = 0.39. Al igual que en los mapas discretos, λ > 0 caracteriza el caos del sistema. Notemos que si el sistema es caotico resulta dificil obtener la soluci´on numericamente, ya que los errores debido a la finitud de ∆t se amplifican exponencialmente.1 Otro metodo de clasificar el caos consiste en generar mapas discretos a partir de la soluci´on de la ecuaci´on diferencial. Estos mapas discretos se analizan mediante las t´ecnicas explicadas en el cap´ıtulo anterior.

7.0.1

Corte de Poincar´ e

Este es un metodo que permite reducir una dimensi´on. Considerese como ejemplo el sistema de Lorenz. Consideremos los puntos de la trayectoria x(t), y(t), z(t) en que x(t) = C y dx/dt > 0, es decir, los puntos en que la trayectoria atraviesa cierto plano en direccion positiva. Esto da una secuencia de instantes tn y los correspondientes valores yn = y(tn ) y zn = z(tn ). Entonces con cada variable se puede hacer un mapa yn+1 = M (yn )

o zn+1 = M (zn )

(7.13)

yn + 1 18 16 14 12 12

14

16

18

yn

Figure 7.3: Corte de Poincare del sistema de Lorenz. Izquierda: imagen de libro. Derecha: calculado por mi, con σ = 10, b = 8/3, r = 28, t ∈ (200, 500).

1 He estudiado que variando ∆t de 0.0001 a 0.01 con el algoritmo NDSolve de Mathematica, el valor de λ aumenta hasta 0.55.

94

CHAPTER 7. CAOS EN ECUACIONES DIFERENCIALES ORDINARIAS z(t)

zn+1

40

45 42.5 40 37.5 35 32.5

35 30 25 20 15 10 202

204

206

208

210

t

32.5 35 37.5 40 42.5 45

zn

Figure 7.4: Izquierda: evoluci´on de la componente z(t) del sistema de Lorenz. Derecha: mapa de los maximos sucesivos de z(t) del sistema de Lorenz. Par´ametros: r = 28, σ = 10, b = 8/3. Mapa de los maximos sucesivos En realidad este corte de Poincar´e, en caso del sistema de Lorenz, no permite visualizar bien, pues tanto el dominio y la imagen del mapa son bidimensionales. Esa es la causa de que el mapa, como funci´on de yn , parezca multivaluado. Una alternativa util es la sucesi´on de m´aximos de la variable z(t). El mapa zn+1 = M (zn ) se muestra en la Fig. 7.4 Con este mapa se puede hacer diagrama de bifurcaciones como funci´on de los paramtros r, σ, b, calcular exponentes de Lyapunov, y todo lo demas que se hace con mapas discretos. Comparando con la recta y = x se ve que siempre |dM/dx| > 1 y no es invertible, lo que implica inmediatamente que el mapa, y por tanto el sistema de ecuaciones diferenciales, son ca´oticos.

7.0.2

Mapa estrobosc´ opico

Este consiste en la funci´on xn = z(nT ), siendo T un valor que determina la escala de la variable t. Esa escala de tiempo es discernible en casos como el del p´endulo forzado. d2 x dx +β + sin x = A cos ωt. 2 dt dt

(7.14)

Este sistema se lleva a la forma can´onica (7.1) introduciendo las variables auxiliares y = dx/dt y z = t. dx = y dt dy = −βy − sin x + A cos ωz dt dz = 1 dt

(7.15) (7.16) (7.17)

La variable z = t se introdujo para que el miembro derecho solo dependa de las variables del espacio de fases. Por tanto las trayectorias en el espacio x, y, z no se cortan. La figuras 7.5

95

y A=1.0

x

2

25

50

A=1.0 75 100

125

150

t 2π

-0.5 1 -1 -1

-2

1

2

x

-1.5 -2

-1 -2.5 -3

-2

Figure 7.5: Trayectoria en el espacio de fases del pendulo forzado. y = dx/dt, siendo x el angulo con la vertical. x(0) = 0, y(0) = 0, β = 0.1, A = 1.0, ω = 1.

y

A=1.2

x 25

2

-5

1

-10

50

A=1.2 75 100

125

150

t 2π

-15 -24

-22

-20

-18

x -20

-1 -25 -2

-30 x 5

10

A=1.2 15

20

25

t 12 π

-5 -10 -15 -20 -25 -30

Figure 7.6: Trayectoria en el espacio de fases del p´endulo forzado. y = dx/dt, siendo x el ´angulo con la vertical. x(0) = 0, y(0) = 0, β = 0.1, A = 1.2, ω = 1.

96

CHAPTER 7. CAOS EN ECUACIONES DIFERENCIALES ORDINARIAS

y 7.6 muestran la proyecci´on del atractor en el plano xy de las trayectorias del pendulo, con condiciones iniciales x = y = z = 0. Para A = 1.0 el atractor difiere poco del oscilador lineal, pero para A = 1.2 el atractor es mucho mas complejo. El periodo T = 2π es determinado por la fuerza perturbadora. El diagrama estrobosc´opico para A = 1.0 muestra que la trayectoria es periodica. Para A = 1.2 muestra que el periodo no es 2π. Al hacer T = 12π se ve claramente que la trayectoria es periodica. N´otese la gran diferencia de amplitud entre los atractores para A = 1.0 y A = 1.2, siendo iguales las condiciones iniciales.