N° d’ordre : /2010 - M /CHRépublique Algérienne Démocratique et Populaire Ministère de l’Enseignement Supérieur et de la Recherche Scientifique Université des Sciences et de la Technologie Houari Boumediene Faculté de Chimie Ecole Doctorale Physique Chimie Théorique Chimie Informatique MEMOIRE Présenté pour l’obtention du diplôme de MAGISTER En Chimie Option : Physique Chimie Théorique Chimie Informatique Par : DEKHIRA Azzeddine Sujet : Soutenu le :…………………….., devant le jury composé de : Etude théorique et simulation des cristaux photoniques et leurs applications en chimie et biochimie Mr. A. AIT KACI Professeur à l’USTHB Président M elle . O. OUAMERALI Professeur à l’USTHB Directrice de thèse M elle . D. HAMMOUTENE Professeur à l’USTHB Examinatrice Mr. O. KRACHNI Professeur à l’UFA S Examinateur Mr. T. ALI ZIANE Maitre de Conférences à l’USTHB Examinateur REMERCIEMENTS Ce travail de mémoire a été effectué au sein de l’école doctorale physique chimie théorique chimie informatique, sous la direction scientifique du Professeur Ourida OUAMERALI, responsable de l’équipe 2 du laboratoire de physico-chimie théorique chimie informatique, à la Faculté de Chimie de l’U.S.T.H.B J’adresse mes profonds remerciements à ma directrice de thèse, professeur Ourida OUAMERALI, qui a toujours montré de l’enthousiasme pour mon travail et pour le sujet nouveau que constituent les cristaux photoniques, pour m'avoir confié ce travail et assurer l’encadrement de cette thèse et Je suis très reconnaissant pour le confiance qu’elle m’a accordée, ses conseils judicieux, sa disponibilité et le soutien constant qu’elle m’a prodigué au cours de l’élaboration de ce travail. Je remercie vivement Monsieur le Professeur M. A. AIT KACI qui m’a fait l’honneur d’accepter la présidence du jury de ce mémoire .Je lui exprime toute ma gratitude pour l’intérêt qu'il a porté à ce travail. Je suis reconnaissant à Melle D. HAMMOUTENE, professeur à l’USTHB de m’avoir honoré de sa présence en étant membre de jury. Je la remercie très respectueusement d'avoir accepté de juger ce travail. Je remercie également, Monsieur O. KRACHNI, professeur à l’université Ferhat Abbes Sétif pour ses conseils et ses encouragements et aussi a bien voulu être membre du jury et examiner ce mémoire. J’exprime ma gratitude envers Monsieur T. ALI ZIANE, maître de conférence à l’USTHB pour ses encouragements et d’avoir accepté de faire partie du jury. Mes remerciements sont adressés également à tous les enseignants de l’ECOLE DOCTORALE Physique Chimie Théorique Chimie Informatique qui ont contribué à ma formation. J’associe à mes remerciements mes camarades de l’école doctorale : S. BOUARAB, A. SADI, N. BENSERADJ, S. REZZOUK et A. BOUROUINA pour l’ambiance chaleureuse de travail et pour nos échanges qui n’ont pas toujours été scientifiques. Je remercie aussi vivement les membres de notre équipe de recherche : S. MOUSSI, M.HADJ BEN ALI, D.KHEFFACHE Y. MOUSSAOUI et M. REKHIS pour leur soutien moral. Pour leur amitié jamais démentie, pour leur soutien moral et leurs encouragements constants, je remercie mes amies : K. BABESSE, N. YAHIAOUI, K. DJILANI, H. BELAID, M. CHEHILI, B. DJAALAB, L. MEHAMELI, A. BESSAS, N. KETTAF, A. MEZIOUD, S. CHEHILI et M. MEZIANE Je suis particulièrement reconnaissant à mes parents, mes sœurs et mes frères qui ne ménagent aucun effort pour me soutenir. Sommaire 1 SOMMAIRE Introduction Générale 6 Bibliographie 10 Chapitre I : Généralités sur les cristaux photoniques 11 I.1 Introduction 12 I.2 Définition 13 I.3 Bref historique 13 I.4 Notion de bande interdite photonique 14 I.5 Caractéristiques des cristaux photoniques 15 I.6 Matériaux BIP à défaut 16 I.7 Classes de cristaux photoniques 17 I.7.1 Cristaux photoniques tridimensionnels 17 I.7.2 Cristaux photoniques Bidimensionnels 20 I.7.3 Cristaux photoniques unidimensionnels 22 I.8 Méthodes d’élaboration 22 I.8.1 Méthodes lithographiques 23 I.8.2 Méthodes holographiques 24 I.8.3 Méthodes d’auto-assemblage 25 I.9 Méthodes d’élaboration 26 I.9.1 Papillons 26 I.9.2 Souris de mer « Aphrodita » 27 I.9.3 Opales naturelles 28 Bibliographie 29 Chapitre II : Etude théorique des cristaux photoniques 31 II.1 Introduction 32 II.2 Equations macroscopiques de Maxwell 32 II.3 Analogie Schrödinger-Maxwell 35 II.3.1 Propriétés des modes harmoniques 37 II.3.2 Loi d'échelle 38 II.3.3 Différences et similarités 39 II.4 Théorème de Bloch 40 Sommaire 2 II.4.1 Réseau direct et réseau réciproque 40 II.4.2 Zones de Brillouin 41 II.4.3 Zone de Brillouin irréductible 41 II.5 Diagramme de bandes 42 II.6 Carte des bandes interdites 44 II.7 Bandes permises et interdites 45 II.7.1 Etude quantique 45 II.7.2 Etude électromagnétique 50 II.8 Vitesse de groupe et vitesse de phase 53 II.9 Conclusion 55 Bibliographie 56 Chapitre III: Méthodes de simulation numérique 57 III.1 Introduction 58 III.2 Méthode de décomposition en ondes planes 59 III.3 Méthode des différences finies dans le domaine 61 III.4 Méthode des éléments finis 64 III.5 Méthode rigoureuse des ondes couplées 66 III.6 Méthode de la ligne de transmission 67 III.7 Méthode des matrices de transfert 69 III.8 Approches hybrides 69 Bibliographie 71 Chapitre IV: Méthode des différences finies dans le domaine temporel 73 IV.1 Introduction 74 IV.2 Equations de Maxwell dans l’espace cartésien 75 IV.3 Réduction à deux dimensions 76 IV.3.1 Polarisation TE 77 IV.3.2 Polarisation TM 77 IV.3.3 Propagation off-plane 77 IV.4 Réduction à une dimension 78 IV.5 Algorithme de Yee 79 IV.5.1 Principe des différences finies centrées 79 IV.5.2 Discrétisation des équations de Maxwell 81 Sommaire 3 IV.5.3 Equations de Maxwell aux différences centrées 85 IV.5.4 Dispersion numérique 87 IV.5.5 Critères de convergence et de stabilité de l’algorithme 89 IV.6 Sources et signaux d'excitation 90 IV.6.1 Impulsion Gaussienne 90 IV.6.2 Excitation sino-gaussienne 91 IV.6.3 Excitation par une onde plane 92 IV.7 Conditions d’absorption aux limites 93 IV.7.1 Bref état de l’art 94 IV.7.2 Conditions périodiques aux limites 95 IV.7.3 Conditions d’Engquist-Majda-Mur 96 IV.7.4 Couches parfaitement adaptées « PML » 98 IV.8 Implémentation des milieux dispersifs 101 IV.8.1 Méthode RC 101 IV.8.2 Méthode ADE 103 IV.9 Conclusion 104 Bibliographie 105 Chapitre V: Méthode de décomposition en ondes planes 107 V.1 Introduction 108 V.2 Equation de Helmholtz 108 V.3 Structure de bandes des cristaux photoniques unidimensionnels 109 V.3.1 Position du problème 109 V.3.2 Calcul de structure de bandes 110 V.3.3 Solution du problème aux valeurs propres 113 V.3.4 Algorithme de la méthode PWE 114 V.4 Structure de bandes des cristaux photoniques 2D et 3D 114 V.4 .1 Cas d'un cristal photonique 3D 114 V.4 .2 Cas d'un cristal photonique bidimensionnel 116 V.5 Développement de Fourier de la fonction diélectrique 117 V.6 Structure de bandes «off-plane » d’un cristal photonique 2D 118 V.7 Structure de bandes d’un cristal photonique avec défaut 118 V.8 Conclusion 121 Sommaire 4 Bibliographie 122 Chapitre VI: Conception et développement d'un logiciel de simulation 123 VI.1 Introduction 124 VI.2 Description et architecture du logiciel 124 VI.2.1 Module d’interface Windows 126 VI.2.2 Module d’entrée 126 VI.2.3 Module de Sortie 126 VI.2.4 Module de Simulation 126 VI.3 Simulateur FDTD 126 VI.4 Simulateur PWE 129 VI.5 Interface graphique 130 VI.6 Validation du module de simulation 132 VI.6.1 Validation du simulateur FDTD 132 VI.6.2 Validation du simulateur PWE 133 VI.7 Conclusion 135 Bibliographie 136 Chapitre VII: Applications des cristaux photoniques 137 VII.1 Introduction 138 VII.2 Capteurs à base de cristaux photoniques 139 VII.2.1 Description de la structure 140 VII.3.2 Modélisation de la structure 140 VII.3 Spectroscopie ultrarapide à base de cristaux photoniques 141 VII.3.1 Lasers à impulsions femtoseconde 141 VII.3.2 Description de la structure 142 Bibliographie 143 Conclusion générale et perspectives 145 Annexes 149 Annexe I : Transformée de Fourier 150 Annexe II : Origine de la bande interdite photonique 152 Glossaire 5 FEM : Finite Elements Method RCWA : Rigorous Coupled Wave Analysis TLM : Transmission Line Matrix Liste des abréviations utilisées CP : Cristal Photonique BIP : Bande Interdite Photonique BPG : Photonic Band Gap 1D : Unidimensionnel 2D : Bidimensionnel 3D : Tridimensionnel TE : Transverse Electrique TM : Transverse Magnétique ZB : Zone de Brillouin PWE : Plane Wave Expansion FDTD : Finite Difference Time Domain PML : Perfectly Matched Layer CPML : Convolution Perfectly Matched Layer UPML : Uniaxial Perfectly Matched Layer TTM : Transfer Matrix Method TF/SF: Total Field / Scattered Field ABC : Absorbing Boundary Conditions RBC : Radiation Boundary Conditions PBC : Periodic Boundary Conditions CPL : Common Photonic Layer RC : Recursive Convolution ADE : Auxiliary Differential Equation GUI : Graphical User Interface FFT : Fast Fourier Transform MEMS : Micro Electro Mechanical Systems Introduction Générale Introduction Générale 7 Ces quinze dernières années, les cristaux photoniques ou matériaux à bande interdite photonique (BIP) [1] ont suscité un intérêt important dans la communauté scientifique. Cet intérêt pour ces matériaux est dû au fait qu’ils ont des propriétés optiques uniques. Les cristaux photoniques sont des matériaux hétérogènes artificiels ou naturels dont l’indice de réfraction varie périodiquement dans les différentes directions de l’espace et constituent à l’heure actuelle une nouvelle classe de matériaux. À l'image des électrons dans les semi-conducteurs, les photons y sont répartis en bandes de transmission séparées par des bandes d'énergies interdites. Cette analogie [2] permet d'envisager l'utilisation des cristaux photoniques pour stocker, localiser, filtrer ou bien guider la lumière. Le développement de ce nouveau type de matériau a ouvert la voie à un nouveau champ de recherche et à des possibilités d’applications très diverses. Cependant, le développement de ces applications se heurte encore à la difficulté rencontrée pour la fabrication et la caractérisation de ces matériaux notamment aux fréquences optiques. Cette difficulté rend coûteuses en temps et argent les études expérimentales systématiques. Il a donc été nécessaire de disposer d’une modélisation théorique et numérique efficace et rapide de ces cristaux permettant d’orienter la fabrication vers des cristaux performants. L’objectif principal de ce travail était de mener une étude théorique et numérique des cristaux photoniques et de développer un logiciel de simulation et d’analyse de ces matériaux, dopé d’une interface graphique et possède une structure modulaire. Cet outil informatique est destiné à remplir la nécessité pour un programme personnalisable et extensible qui peut satisfaire aux besoins spécifiques de la recherche dans le domaine de photoniques. Dans le premier chapitre, des notions et des concepts de base sur les cristaux photoniques sont présentés. Pour cela, la description et l’historique de matériaux à bande interdite photonique, la notion de bande interdite photonique et les différents types de défauts intentionnels sont exposés. Ensuite un aperçu de différentes classes de ces matériaux ainsi que leurs propriétés physiques est donné. Nous continuerons par une présentation des méthodes de fabrication des cristaux photoniques artificiels et une description de quelques exemples de matériaux BIP naturels. Introduction Générale 8 Le deuxième chapitre est consacré à l’étude théorique des cristaux photoniques. L'idée principale consiste à exploiter l’analogie entre les semi-conducteurs électroniques, dont la périodicité atomique interdit la propagation des électrons dans certaines bandes d’énergie, et les photons piégés dans des structures photoniques et par conséquent, l’analogie entre les équations de Maxwell sous leurs formes fréquentielles et celle de Schrödinger. L’équation de propagation obtenue, qui représente un problème aux valeurs propres, est résolue en utilisant des concepts et des outils développés en physique du solide et en mécanique quantique tels que le théorème de Bloch [3, 4] et la transformée de Fourier. Le troisième chapitre est consacré aux rappels de quelques méthodes qui peuvent être mises en œuvre pour la modélisation des cristaux photoniques. Nous mettons l'accent sur les méthodes les plus utilisées telles que la FDTD [5, 6], La méthode PWE [7], la méthode des éléments finis [8] et celle de la ligne de transmission [9] ainsi que les méthodes hybrides. Le quatrième chapitre présente la formulation et l’implémentation informatique de la méthode de différences finis dans le domaine temporel « FDTD » basé sur la discrétisation des équations de Maxwell exprimées en coordonnées cartésiennes, pour la modélisation des matériaux à bande interdite photonique. Cette méthode permet essentiellement de simuler la propagation de la lumière dans les structures photoniques finies et d’obtenir les coefficients de réflexion et de transmission. Les conditions d’absorption aux limites « PML » [10] pour décrire l’espace libre sont intégrées aux codes FDTD ainsi que les deux modèles de Debye et de Lorentz pour décrire les milieux dispersifs. Le cinquième chapitre présente la formulation et l’implémentation de la méthode de décomposition en ondes planes « PWE » qui consiste à résoudre l’équation d’onde dans l’espace fréquentiel en développant le champ magnétique sur une base d’ondes planes. Cette méthode permet de calculer les diagrammes de bandes de cristaux photoniques. Pour traiter des structures photoniques en présence de défauts, la technique de supercellules [11] est intégrée aux codes PWE. Le sixième chapitre montre d’une manière générale les organigrammes des algorithmes développés ainsi que les outils utilisés pour la réalisation d’un logiciel de simulation et d’analyse des cristaux photoniques avec une interface graphique (GUI). Ce logiciel est développé en langage orienté objet C++Builder avec une structure modulaire, personnalisable et extensible. Le cœur du logiciel est le module de simulation qui consiste en deux simulateurs, l’un est développé à la base de la méthode FDTD et l’autre à la base de la Introduction Générale 9 méthode PWE. Plusieurs tests de validation sur des structures photoniques en niobate de lithium « LiNbO3 », Arséniure de Galium « GaAs », Dioxyde de titane « TiO2 » et en Silicium sont effectués pour évaluer la fiabilité du logiciel. Le septième chapitre est consacré aux applications des cristaux photoniques. Deux applications d’intérêt chimique et biochimique sont décrites : capteurs chimiques et biochimiques à base de cristaux photoniques [12] et Spectroscopie femtoseconde [13]. Introduction Générale 10 Bibliographie [1] E. Yablonovitch, Inhibited spontaneous emission in solid-state physics and electronics. Phys. Rev. Lett. 58, 2059–2062 (1987). [2] G. Malpuech, A. Kavokin, G. Panzarini, and A. Di Carlo, Theory of photon Bloch oscillations in photonic crystals, Physical Review B 63, 035108 (2001). [3] C. Kittel, Introduction to Solid State Physics, John Wiley & Sons, Inc. 1996. [4] R. D. Meade, A.Devenyi, J. D. Joannopoulos, O. L. Alerhand, D. A. Smith et K. Kash, Novel applications of photonic band gap materials: Low loss bends and Q cavities, [5] K. S. Yee, Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media, IEEE Transactions on Antennas and Propagation, vol. 14, no. 3, pp. 302{307, 1966. [6] A.Taflove and M.E.Brodwin, IEEE Transactions on Microwave theory and Technique, MTT-23, No 8, August 1975. [7] K. M. Ho, C. T. Chan, and C. M. Soukoulis, Existence of a photonic gap in periodic dielectric structures, Phys. Rev. Lett., vol. 65, no. 25, pp. 31523155, 1990. [8] R. L. Courant, Variational Methods for the Solution of Problems of Equilibrium and Vibration, Bulletin of the American Mathematical Society 49: 1-23., 1943. [9] P. B. Johns et R. L. Beurle, Numerical solution of 2-dimensional scattering problems using a transmission-line matrix, Proceedings IEE, vol. 118, p. 1203–1208, sept. 1971. [10] J.P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics 114, p. 185 (1994). [11] A. Yariv, Y. Xu, R. K. Lee, and A. Scherer Coupled-resonator optical waveguide: a proposal and analysis, Optics Lett 24, 711 (1999). [12] T. Stomeo, M. Grande, A. Qualtieri, A.Passaseo, A.Salhi, M. Vittorio. Microelectronic Engineering, 2007, Vol. 84, issue 5-8, pp 1450-1453 (2007). [13] C. Lecaplain, A. Hideur, S. Février, P. Roy, Mode-locked Yb-doped Bragg fiber laser, Optics Letters, Vol. 34, no. 18, pp.2879-2881 (2009). CHAPITRE I Généralités sur les cristaux photoniques CHAPITRE I : Généralités sur les cristaux photoniques 12 I.1 Introduction Depuis une décennie, une communauté de chercheurs rassemblant opticiens, physiciens et chimistes s’est fixé l’objectif ambitieux de réaliser un matériau qui serait, pour les photons, l’analogue de ce qu’est un cristal semi-conducteur pour les électrons. Cette nouvelle classe de matériaux a suscité un très vif intérêt dans le monde de la recherche et ceci dans plusieurs secteurs de la physique et de la chimie. Il s'agit des structures périodiques diélectriques ou métalliques, rencontrées sous les appellations «cristaux photoniques» ou « matériaux à bande interdite photonique », qui présentent des états photoniques structurés en bandes interdites et passantes. En effet, dans un cristal semi-conducteur la périodicité atomique empêche les électrons de prendre n’importe quelle valeur d’énergie ; elle doit appartenir à certaines gammes d’énergies séparées par des « bandes d’énergies interdites ». Ces sont encore appelées bandes interdites électroniques (electronic band gap). Tout l’intérêt des semi- conducteurs découle de l’existence de cette zone. E. Yablonovitch démontra [1], dans le but de contrôler directement l’émission de lumière, la possibilité théorique de fabriquer dans des matériaux diélectriques, des structures qui possèdent une périodicité semblable à celles des cristaux atomiques. Dans ces structures, les bandes interdites ne concerneraient plus les électrons mais les photons. C’est ainsi que naquirent les concepts de Bande Interdite Photonique (BIP), en anglais Photonic Band Gap (PBG) et de cristaux photoniques (CP). Les cristaux photoniques existent dans la nature à l’état minéral et biologique. Les opales sont des minéraux composés d’arrangements de sphères de silice hydratée et l’origine de la coloration de nombreuses espèces animales et végétales provient aussi de motifs périodiques. Toutefois, les cristaux photoniques qui constituent un domaine de recherche très dynamique sont souvent le résultat de synthèses artificielles. Les progrès récents dans les techniques et les méthodes de fabrication des cristaux photoniques permettent de réaliser des structures à l’échelle du nanomètre qui contrôlent la lumière visible et infrarouge et permettent ainsi d’envisager de nombreuses applications potentielles révolutionnaires. CHAPITRE I : Généralités sur les cristaux photoniques 13 I.2 Définition Les cristaux photoniques, encore appelés matériaux à bandes interdites photoniques (BIP), sont des matériaux diélectriques, semi-conducteurs ou métalliques artificiellement ou naturellement structurés dont la constante diélectrique varie périodiquement à l’échelle de la longueur d’onde selon une ou plusieurs directions de l’espace. Par analogie avec la bande d’énergie interdite électronique caractérisant les réseaux cristallins atomiques, les structures photoniques possèdent une bande de fréquences interdites dans laquelle aucune onde électromagnétique ne peut se propager, indépendamment de la polarisation et de la direction de propagation. Cette propriété intéressante offre aux cristaux photoniques la possibilité du contrôle de la propagation sans absorption des ondes électromagnétiques et permet ainsi des perspectives nouvelles pour la manipulation de la lumière. I.3 Bref historique Malgré le fait que ce n'est que pendant les dernières décennies que les cristaux photoniques ont attiré une grande attention, les premières hypothèses sur la possibilité de contrôler la propagation de la lumière utilisant des structures périodiques se rapportent à 1887 avec les travaux de Lord Rayleigh [2]. En 1972, une étude théorique détaillée de structures optiques unidimensionnelles a été réalisée par V.P. Bykov [3]. Il a été le premier à examiner l'effet de bandes interdites sur l'émission spontanée provenant d'atomes et de molécules intégrées à la structure. Bykov fit aussi des hypothèses sur l'emploi de structures bidimensionnelles et tridimensionnelles. On considère souvent que le domaine des cristaux photoniques à démarré en 1987, quand E. Yablonovitch et S. John [4] ont introduit, séparément et dans des contextes différents, le concept de matériaux à bandes interdites photoniques. La motivation principale de Yablonovitch était d'appréhender la densité d'états photoniques, par analogie à la densité d'états électroniques, dans le but de contrôler l'émission spontanée de matériaux intégrés aux cristaux photoniques. John, quant à lui, voulait utiliser les cristaux photoniques pour modifier la localisation et le contrôle de la lumière. CHAPITRE I : Généralités sur les cristaux photoniques 14 Après 1987, le nombre de publications concernant les cristaux photoniques commença à croître exponentiellement. Cependant, à cause de la difficulté de fabrication de ces structures pour qu'elles soient effectives dans le spectre visible, les premières études étaient soit théoriques, soit dans les micro-ondes. En 1991, A. Genack et al [5] ont montré expérimentalement l’existence de l’effet de localisation de la lumière dans les structures périodiques. En même temps, Yablonovitch et al [6] ont démontré expérimentalement la possibilité de réaliser une structure diélectrique capable de réfléchir la totalité d’un rayonnement électromagnétique, quelle que soit la direction incidente et dans le domaine des micro-ondes. En 1993, Yablonovitch conçoit le premier cristal photonique tridimensionnel possédant une bande interdite dans les micro-ondes. Ce cristal photonique s’appelle d’après son inventeur « la Yablonovite ». En 1996, Thomas Krauss fit la première démonstration d'un cristal photonique bidimensionnel dans le spectre du visible [7]. Cela ouvrit la voie à la fabrication de cristaux photoniques par les méthodes utilisées dans le secteur des semi-conducteurs. En 1998, l'opale inverse artificielle a été obtenue expérimentalement [8]. Le diamètre de la sphère dans la structure était d'environ 1 µm, et la distance entre les sphères est très faible. L'indice de réfraction du matériau utilisé ( TiO 2 ) entre les sphères est de 2.8. En 2000, le premier cristal photonique tridimensionnel avec une bande interdite photonique complète dans le domaine infrarouge proche a été obtenu [9]. Un tel cristal photonique s'est composé des sphères de silicium disposées dans une maille d’un cristal de type « diamant ». Au cours des dernières années, la recherche dans le domaine des cristaux photoniques a connu une expansion extraordinaire et a couvert presque toutes les disciplines scientifiques en réalisant des progrès sans précédent. I.4 Notion de bande interdite photonique Dans un semi-conducteur, la variation périodique du potentiel d’interaction entre électrons et atomes fait que les électrons n’ont accès qu’à certains niveaux d’énergie, des bandes d’énergie permises, séparées entre elles par des bandes d’énergies interdites. Ce concept de bandes permises et interdites peut être étendu au comportement des photons dans CHAPITRE I : Généralités sur les cristaux photoniques 15 un cristal photonique. A cause de la variation périodique de l’indice de réfraction dans un cristal photonique, l’énergie des photons est quantifiée en bandes permises et en bandes interdites, appelées aussi gaps. Les bandes permises et interdites d’un CP se regroupent dans un diagramme de bandes photoniques, qui est une représentation des fréquences possibles pour l’onde électromagnétique au sein du CP en fonction de son vecteur d’onde. En revanche, dans le cas d’un CP présentant une symétrie cristalline adaptée, un contraste d’indice de réfraction suffisamment élevé et constitué de motifs élémentaires de forme appropriée ; les bandes interdites peuvent devenir assez larges et se recouvrir pour une certaine gamme de fréquences. La propagation de la lumière est de la sorte interdite dans le matériau pour ces fréquences, selon toutes les directions de l’espace. On parle alors de « bande interdite photonique complète ». Une bande interdite photonique d'un cristal est dite complète (ou totale) lorsque, pour le domaine de fréquences considéré, le cristal ne supporte aucun mode électromagnétique de propagation ; c'est-à-dire qu'une onde dont la fréquence est dans la bande interdite totale ne peut pas se propager dans le cristal quelques soient sa polarisation et sa direction de propagation. Expérimentalement, une bande interdite est mise en évidence en mesurant la réponse du matériau soumis à un faisceau lumineux, en transmission ou en réflexion. Cette dernière est caractérisée par l’apparition d’un minimum de la transmission et par conséquent un maximum de la réflexion. I.5 Caractéristiques des cristaux photoniques La dimensionnalité Elle est déterminée par la périodicité de l’indice de réfraction. La périodicité d’un cristal photonique peut s’étendre à une, deux ou trois dimensions. La symétrie La position des éléments d’un CP détermine la symétrie du réseau. Par exemple, pour un CP 3D de particules sphériques, une symétrie cubique, hexagonale compacte (hc) ou cubique à face centrée (cfc) peut être obtenue. CHAPITRE I : Généralités sur les cristaux photoniques 16 La topologie La topologie rend compte de l’architecture, de la compacité du matériau. Un réseau d’une symétrie donnée peut présenter des topologies différentes (cas de briques constitutives interpénétrées, en contact ou isolées). Le paramètre du réseau C’est la distance fondamentale entre deux éléments constitutifs. Il détermine la région spectrale où le CP interagit avec l’onde électromagnétique. Le contraste d’indice de réfraction Ce paramètre est défini comme le rapport n1/n2 entre les indices de réfraction des éléments et de la matrice. Il offre une idée générale de la force de diffusion des deux matériaux composants du cristal photonique. I.6 Matériaux BIP à défaut Toujours par analogie avec les cristaux semi-conducteurs, les fonctionnalités des cristaux photoniques peuvent être exaltées en insérant volontairement et de façon contrôlée des défauts au sein de leur structure. On parle alors de défauts extrinsèques, en opposition avec des défauts non intentionnels, intrinsèques aux CP, comme des imperfections dans la structure apparaissant lors de leur fabrication. Ces derniers étant présents dans les matériaux de façon aléatoire, ils entraînent une dégradation des propriétés optiques et sont donc nuisibles aux applications finales. La création d’un défaut extrinsèque est causée par la rupture de la périodicité de l’indice de réfraction. Comme pour les semi-conducteurs, où des niveaux d’énergie apparaissent dans le gap lors de l’insertion d’impuretés (atomes autres que ceux du cristal), les défauts extrinsèques au sein de CP créent des niveaux d’énergies permis, nommés « modes de défauts », pour des fréquences particulières dans la bande interdite. Il existe deux principaux types de défauts : les défauts ponctuels et les défauts étendus. Les premiers, associés à une rupture locale de périodicité, se traduisent par la présence de modes électromagnétiques à des fréquences discrètes, analogues aux défauts électroniques. Les seconds, que l’on peut considérer comme analogues aux dislocations, peuvent donner lieu à des bandes permises de propagation, là où se trouve une bande interdite dans le cristal idéal. CHAPITRE I : Généralités sur les cristaux photoniques 17 L'insertion de défauts dans les cristaux photoniques nécessite une modification contrôlée d’un ou plusieurs paramètres au cours du processus de fabrication. Les paramètres les plus considérés pour cette opération sont les suivants: Dimensions des motifs élémentaires Pour rompre la périodicité d’une structure BIP, on peut modifier la taille du motif élémentaire (Fig. 1.1) qui compose le cristal photonique. Figure 1.1 : Défaut de dimension du motif élémentaire Distance entre motifs élémentaires On peut aussi jouer sur l’espace qui existe entre les motifs élémentaires des réseaux cristallins (Fig. 1.2). Figure 1.2 : Défaut de distance entre motifs élémentaires Valeur de la permittivité relative des motifs élémentaires Il est possible de changer localement la nature du matériau et plus concrètement, changer la valeur de la permittivité relative (Fig. 1.3). Figure 1-3 : Défaut sur la permittivité relative CHAPITRE I : Généralités sur les cristaux photoniques 18 Défaut par vacuité Le défaut par vacuité correspond à l’élimination de motifs élémentaires qui se trouvent remplacés par la permittivité de fond (Fig. 1.4). Figure 1.4 : Défaut par vacuité I.7 Classes de cristaux photoniques Il existe différents types de cristaux photoniques, à classer selon leur dimensionnalité. A une dimension (1D), on retrouve les bien connus miroirs de Bragg (Fig. 1.5a) formés d'une alternance de couches de bas et haut indice. Le principe des miroirs de Bragg peut être généralisé à 2 ou 3 dimensions, constituant des cristaux photoniques bidimensionnels (Fig. 1.5b) ou tridimensionnels (Fig. 1.5c). I.7.1 Cristaux photoniques tridimensionnels Les cristaux photoniques 3D ont attiré et attirent toujours de nombreux efforts de recherche. Ils constituent la seule structure qui permet d'obtenir une bande d'énergie interdite dans toutes les directions de l'espace. Le premier cristal photonique 3D a été fabriqué par K.M Ho et al. [10]. Il était formé de sphères de silicium arrangées sur une structure diamant. Mais Figure 1.5: Schéma de cristaux photoniques 1D, 2D ou 3D. Les différentes couleurs représentent des matériaux de constants diélectriques différents. (1D) (2D) (3D) CHAPITRE I : Généralités sur les cristaux photoniques 19 l'histoire retient généralement la célèbre Yablonovite, structure 3D pour les micro-ondes fabriquée en 1993 par E. Yablonovitch. Au fil des années, les scientifiques ont cherché à réduire la dimension des motifs, en utilisant plusieurs méthodes, pour aboutir aujourd’hui à des cristaux photoniques présentant une bande interdite dans le proche infrarouge et le visible. De nombreuses structures tridimensionnelles ont été proposées. Les deux suivantes ont attiré le plus d'efforts de recherche: Structures « Tas de bois » : Ces structures 3D sont obtenues en déposant par couches successives des rubans de silicium polycristallin dans des tranchées de silice. Après avoir bâti la structure, la silice est retirée pour obtenir un cristal photonique 3D Si/air dont le contraste d'indice est suffisant pour ouvrir une bande d'énergies interdites omnidirectionnelle [11] (Fig. 1.6). Des cristaux photoniques semblables ont été fabriqués sur GaAs par Noda et al. [12] par un procédé de fusion/élimination du substrat. Cette technique utilise des technologies standards de micro-fabrication des semi-conducteurs et permet l'introduction déterministe de défauts dans les cristaux fabriqués. Figure 1.6 : Image MEB (Microscopie Electronique à Balayage) d’un cristal photonique du type tas de bois Opales : Ces structures forment une famille originale de cristaux photoniques 3D. Elles sont obtenues chimiquement par auto-assemblage (Fig. 1.7). La première opale a été obtenue par sédimentation de sphères de silice (SiO2) en solution: la gravité arrange ces sphères selon un réseau cubique à faces centrées [13]. CHAPITRE I : Généralités sur les cristaux photoniques 20 Le nombre important de défauts dans les premières opales a été fortement réduit grâce à des techniques de croissance auto-organisées proposées par Y.A. Vlasov [14]. La plupart de ces cristaux colloïdaux ne présentent pas de bandes d’énergie interdites, à cause du faible contraste d’indice. Cependant, ces structures servent d’empreinte pour la réalisation d’opales inverses à partir de l’infiltration d’un matériau de haut indice (Fig. 1.8). Les sphères initiales sont ensuite dissoutes pour aboutir à structure finale de sphères d’air dans une matrice de haut indice. Figure 1.7 : Vue de MEB d’une opale artificielle directe Figure 1.8 : Vue de MEB d’une opale artificielle inverse I.7.2 Cristaux photoniques Bidimensionnels Les difficultés de fabrication des structures 3D ont conduit à envisager la réalisation et l'étude de structures 2D. Un cristal photonique 2D parfait est périodique dans le plan (Oxy) et infiniment long dans la direction (Oz). Il possède une bande interdite dans le plan (Oxy). Ces systèmes n'existent pas dans la réalité mais de bonnes approximations peuvent être obtenues. CHAPITRE I : Généralités sur les cristaux photoniques 21 Figure 1.9 : Exemple de cristal photonique bidimensionnel imagé en microscopie électronique L'insertion de défauts est plus simple que dans les cristaux photoniques 3D. Pour compenser l'absence de bande interdite dans la direction perpendiculaire au plan de périodicité des cristaux 2D, la lumière peut être confinée dans une hétérostructure d'indice. Cette dernière se compose généralement d'une couche de diélectrique entourée de deux autres couches diélectriques d'indices de réfraction plus faibles. A deux dimensions, il est nécessaire de considérer deux polarisations différentes: TE (avec le champ E perpendiculaire à l'axe des trous) et TM (où E est parallèle à l'axe des trous). Ces deux polarisations sont découplées et donnent lieu à deux diagrammes de bande indépendants. Il n'existe donc pas forcément une bande interdite dans les deux cas. Il existe de nombreux degrés de liberté lors de la conception d'un cristal photonique 2D. En particulier, il est possible pour un type de réseau choisi d'ajuster le paramètre de maille et le facteur de remplissage surfacique (rapport air/surface totale). Ces paramètres influencent directement les propriétés et l'allure du diagramme de bandes associé au cristal photonique réalisé, en particulier la largeur et la position de la bande interdite. La configuration la plus propice à l'obtention d'une bande interdite complète (c'est-à-dire en TE et en TM) est le réseau triangulaire de trous dans un diélectrique de haut indice de réfraction. CHAPITRE I : Généralités sur les cristaux photoniques 22 I.7.3 Cristaux photoniques unidimensionnels Les cristaux photoniques 1D sont les plus simples à réaliser. Ils sont obtenus en empilant périodiquement des couches planes de diélectriques d’indices de réfraction différents. À chaque interface entre deux couches, la lumière est partiellement réfléchie et transmise. Selon la valeur des déphasages (qui eux-mêmes dépendent de la longueur d’onde) on obtient des interférences destructives ou constructives. Les interférences constructives des ondes réfléchies entraînent une réflexion totale. Ainsi, pour certaines longueurs d’onde, la structure multicouche se comporte comme un miroir. Figure 1.10 : Cristal photonique unidimensionnel (miroir de Bragg) Cette réflectivité est la manifestation d’une bande interdite photonique. Cependant, dans les cristaux photoniques 1D la lumière monochromatique n’est réfléchie que lorsqu’elle se propage dans une direction proche de la normale à la structure multicouche. I.8 Méthodes d’élaboration La théorie derrière les propriétés optiques des cristaux photoniques a été largement étudiée au cours des deux dernières décennies et un certain nombre de phénomènes fascinants ont été prévus. Mais la réalisation expérimentale des structures nécessaires pour tester ces prédictions a fait défaut dans de nombreux cas. La raison en devient évidente si l'on prend en compte le fait que la plage de travail pour un cristal photonique est dictée par la périodicité spatiale de son indice de réfraction. Par conséquent, si l'on veut opérer dans la partie visible CHAPITRE I : Généralités sur les cristaux photoniques 23 ou proche infrarouge du spectre électromagnétique, des modulations spatiales de l'indice de réfraction de quelques centaines de nanomètres à un micron sont nécessaires. Cela représente un défi considérable pour la technologie actuellement disponible. En principe, on aurait envie d’avoir une technique efficace qui soit facile à mettre en œuvre, faible coût, et qui conduise à des structures reproductibles de bonne qualité impliquant un délai raisonnable. En ce sens, un certain nombre de méthodes de fabrication ont été inspirées ou directement empruntés à d'autres disciplines. La plupart de ces méthodes peuvent être divisés en trois groupes, chacun ayant des avantages et des inconvénients : méthodes lithographiques, holographiques et d’auto-assemblage. I.8.1 Méthodes lithographiques. Les techniques lithographiques sont fréquemment utilisées dans le domaine de la microélectronique pour fabriquer des composants électroniques. Une couche de résine photo ou électro-sensible telle que le PMMA (Polyméthylmétacrylate) est déposée sur un matériau d’indice de réfraction élevé. La procédure débute par l’enregistrement d’un réseau 2D dans le matériau suivant un procédé de photolithographie ou de lithographie électronique. Dans les zones irradiées, la résine est fragilisée par le rayonnement et éliminée par un solvant ; alors que dans les zones non traitées, elle demeure intacte et protège le substrat. Le motif dessiné dans la résine est transféré dans la matière, par une étape de gravure. On utilise généralement des ions qui viennent frapper la surface et creusent la matière jusqu'à une profondeur voulue. Afin de former un CP 3D, des couches de semi-conducteurs ainsi structurées sont empilées les unes sur les autres. Figure 1.11 : Clichés de MEB de CP élaborés par photolithographie CHAPITRE I : Généralités sur les cristaux photoniques 24 La figure 1.11 montre un exemple d’un cristal photonique élaboré par photolithographie [15]. Dans cette structure, dite en « tas de bois », chaque couche est formée par des bâtonnets parallèles et orientés à 90° par rapport à ceux de la couche sous-jacente, de sorte que les points de contact forment une structure diamant. Cette méthode offre aussi la possibilité de réaliser des défauts intentionnels. Les méthodes lithographiques permettent l’élaboration de CP 3D d’architecture hautement contrôlée, mais limitée à quelques couches. Le procédé est particulièrement onéreux et les nombreuses étapes de fabrication nécessitent un temps considérable. I.8.2 Méthodes holographiques. Le principe de l’holographie [16] consiste à enregistrer l’hologramme créé par l’interférence entre plusieurs faisceaux lumineux cohérents dans une résine photosensible. Afin de créer une structure 3D, quatre sources lumineuses sont requises. La partie de la résine non exposée est dissoute, révélant une structure 3D dont la périodicité et la symétrie sont parfaitement contrôlées par des paramètres expérimentaux comme l’intensité des lasers. La figure 1.12 montre un exemple de CP 3D élaboré suivant cette méthode avec la résine photosensible SU8.19 Cette méthode présente de nombreux avantages. Le temps d’élaboration est très court (quelques minutes), de nombreuses symétries sont accessibles, le procédé est relativement bon marché et adapté pour une production à grande échelle. Enfin, l’addition de défauts optiquement actifs dans la résine avant l’exposition aux lasers est réalisable [17]. Figure 1.12 : Clichés de MEB de CP élaborés par holographie CHAPITRE I : Généralités sur les cristaux photoniques 25 I.8.3 Méthodes d’auto-assemblage L'approche la plus populaire pour la fabrication des cristaux photoniques 3D est celui de l'auto-assemblage [18]. Cette méthode est basée sur la tendance naturelle des particules colloïdales monodisperses de s'auto-assembler dans des rangées organisées communément appelées opales artificielles. Son potentiel en tant que cristaux photoniques a déjà été reconnu dans la première proposition en 1987 et peu après les premières opales artificielles ont été caractérisées optiquement en termes de bandes photoniques. On utilise alors un très grand nombre de particules identiques en interaction et on laisse le système former spontanément une structure qui s’organise de l’échelle microscopique à l’échelle macroscopique. Pour former des cristaux photoniques, on utilise donc généralement cette auto-organisation vers un état méta-cristallin à l’équilibre thermodynamique d'une suspension de particules. On cherche à influer sur ces assemblages en contrôlant des paramètres macroscopiques telles que la température, la concentration en particules, leurs charges de surface, la viscosité du solvant de la suspension. Il existe de nombreuses méthodes permettant de construire des structures par auto-assemblage à 2D et 3D telles que la sédimentation et le dépôt de Langmuir-Blodget… Cependant l'auto-assemblage ne peut conduire qu'aux structures qui sont thermodynamiquement les plus stables, c'est-à-dire les structures compactes, hexagonales et cubiques à faces centrées ou encore à un mélange aléatoire des deux. D’un point de vue optique, ces structures à fortes compacités et hautes symétries présentent peu d’intérêt car par exemple, pour l’ouverture d’une bande interdite photonique; les structures cubiques faces centrées et hexagonales compactes étaient défavorables avec des billes de silice. Les échantillons initiaux construits par la sédimentation, a eu un certain nombre d'inconvénients quant à la difficulté de contrôler l'épaisseur de l'échantillon et le fait qu'ils ne sont pas faciles à manipuler. Ces inconvénients ont ensuite été éliminés par l'introduction de la méthode de dépôt vertical. Parmi les autres inconvénients, on peut citer le contraste faible d’indice de réfraction et la symétrie fixe des échantillons. Une solution pour augmenter le contraste d'indice de réfraction d’une opale artificielle est d'infiltrer ses pores avec un matériau d'indice de réfraction élevé, puis d'enlever la structure originale, pour obtenir ce qu'on appelle une «opale inverse » (Fig. 1.13b). CHAPITRE I : Généralités sur les cristaux photoniques 26 Figure 1.13 : (a) : Les sphères de silicium sont assemblées directement sur le wafer de Si pour former l’opale, (b) : La structure opale est infiltrée avec du silicium puis les sphères de SiO2 sont enlevées par gravure mouillée (opale inverse). I.9 Matériaux à bandes interdites photoniques naturels La nature dispose de nombreux moyens pour produire des effets optiques impressionnants. Il est intuitif d’attribuer les couleurs du monde animal, végétal, minéral à l’absorption sélective de la lumière due à la présence de pigments. Ainsi, en absorbant la lumière rouge et bleue, la chlorophylle donne leur couleur verte aux végétaux. En revanche, certaines couleurs ne peuvent pas être expliquées simplement par un phénomène d’absorption de la lumière. Il existe des structures naturelles qui peuvent avoir des propriétés ayant les mêmes caractéristiques que les cristaux photoniques artificiels. En effet, les colorations vives de certaines espèces sont parfois dues à la présence de structures très complexes, à caractère périodique. Nous allons présenter brièvement certains de ces matériaux. I.9.1 Papillons Les papillons comptent parmi les insectes les plus colorés que nous retrouvons dans la nature. Ils sont en fait très largement tributaires de la lumière et équipés d’un arsenal impressionnant pour gérer cette interaction avec les ondes électromagnétiques [19]. Si on effectue une analyse microscopique des ailes, on découvert que celles-ci sont constituées par des écailles qui ont des structures géométriques dans lesquelles un des paramètres varie en continu. Sur la figure 1.14, apparaît un papillon et la coupe d’une aile. Sur (a) (b) CHAPITRE I : Généralités sur les cristaux photoniques 27 cette coupe observée au microscope électronique, apparaît une structure qui a un comportement de matériau à bande photonique interdite. C’est elle qui donne à certains papillons des couleurs iridescentes. En effet, ce réseau réfléchit la lumière pour certaines longueurs d’onde dans des directions différentes en fonction de la longueur d’onde. Figure 1.14 : La figure à droite présente l’agrandissement d’une aile de papillon. On voit un arrangement périodique des écailles. I.9.2 Souris de mer « Aphrodita » Des scientifiques australiens et britanniques des universités de Sydney et d'Oxford ont trouvé un ver marin possédant des épines qui constituent des cristaux photoniques plus efficaces que ceux fabriqués par l'homme jusqu'à présent [20]. Cet animal au nom charmant : Aphrodita est appelé « souris de mer » de l'anglais « seamouse » (Fig. 1.15a). Cet animal est partiellement recouvert d’épines irisées, elles-mêmes constituées par un arrangement périodique de cylindres creux (Fig. 1.15c). Chaque cylindre ayant un diamètre de l’ordre de la longueur d’onde de la lumière visible. Cette dernière est diffractée par le réseau organisé de cylindres. Figure 1.15 : Photographies d’une souris de mer (a) et d’une de ses épines (b). Coupe d’une épine observée en microscopie électronique à transmission (c). (a) (b) (a) (c) CHAPITRE I : Généralités sur les cristaux photoniques 28 I.9.3 Opales naturelles L’opale est un minéral typique d’origine sédimentaire. Elle se forme par dépôt chimique d’eaux très riches en silice et par accumulation de squelettes d’organismes marins. En effet, c’est un minéral colloïdal amorphe, ou micro cristallin. On la trouve en globules et en croûtes de coloris variés toujours magnifiquement iridescents. Elle contient de la silice et de l’oxygène, dans un rapport de un à deux (comme le quartz) ainsi que de l’eau. Figure 1.16 : Image au microscope électronique d’une opale naturelle constituée d’un réseau quasi-périodique de billes de silice Son étude au microscope électronique a permis de mettre en évidence sa structure. Elle est formée de petites sphères de silice environnées d’espaces vides, équidistants entre eux. Face aux ultraviolets, elle a souvent une fluorescence jaune ou verte. De même, les microbilles de silice peuvent être considérées comme un réseau de diffraction de la lumière incidente. Même si le contraste d’indice entre l’air et la silice (n = 1.5) est faible, on peut utiliser ce genre de structures avec un contraste d’indice important pour réaliser des structures à bande interdite photonique [21] CHAPITRE I : Généralités sur les cristaux photoniques 29 Bibliographie [1] E. Yablonovitch, Inhibited spontaneous emission in solid-state physics and electronics. Phys. Rev. Lett. 58, 2059–2062 (1987). [2] L. Rayleigh, On the maintenance of vibrations by forces of double frequency, and on the propagation of waves through a medium endowed with a periodic structure. Philosophical Magazine 24, 145–159 (1887). [3] V.P. Bykov, Spontaneous emission in a periodic structure. J. Exp. Theor. Phys. 35, 269 (1972). [4] S. John, Strong localization of photons in certain disordered dielectric superlattices. Phys. Rev. Lett. 58, 2486–2489 (1987). [5] A. Genack and N. Garcia, Observation of Photon localisation in a Three-Dimensional Disordered System, Phys. Rev. Lett., 66 (16), pp : 2064, (1991). [6] E. Yablonivitch, T. J. Gmitter, K. M. Leung, Photonic Band Structure: The Face- centred-Cubic Case Employing Nonspherical Atoms, Phy. Rev. Lett., 67, p2295-2298 (1991). [7] T.F. Krauss, R.M. De la Rue et S. Brand, Two dimensional photonic band gap structures operating at near-infrared wavelengths, Nature 383, pp. 699-702, (1996). [8] J.E.G.J. Wijnhoven, L.V. Willem, Preparation of photonic crystals made of air spheres in Titania. Science 281, 802–804 (1998). [9] A. Blanco, E.Chomski, S.Grabtchak, et al. , Large scale synthesis of a silicon photonic crystal with a complete three-dimensional bandgap near 1.5 micrometres. Nature 405, 437–440 (2000). [10] K. M. Ho, C. T. Chan, and C. M. Soukoulis, Existence of a photonic gap in periodic dielectric structures, Phys. Rev. Lett. 65, 25 (1990). [11] S. Y. Lin, J. G. Fleming, D.L. Hetherington, B.K. Smith, R. Biswas, K. M. Ho, M. M. Sigalas, W.Zubrzycki, S.R. Kurtz, and J. Bur, A three-dimensional photonic crystal operating at infrared wavelengths, Nature 394, 6690 (1998). [12] S. Noda, K. Tomoda, N. Yamamoto, and A. Chutinan, Full Three Dimensional Photonic Bandgap Crystals at Near-Infrared Wavelengths, Science 289, 5479 (2000). [13] J. D. Joannopoulos, « Self-assembly lights up » , Nature, vol. 414, no. 6861, pp. 257- 258, 2001. CHAPITRE I : Généralités sur les cristaux photoniques 30 [14] Y.A. Vlasov, X.-Z. Bo, J.C. Sturm, and D.J. Norris, « On-chip natural assembly of silicon photonic bandgap crystals», Nature, vol. 414, no. 6861, pp. 289-293, 2001. [15] Lin, S. Y.; Fleming, J. G.; Hetherington, G. L.; Smith, B. K.; Biswas, R.; Ho, K. M.; Siglas, M. M.; Zubrycki, W.; Kurtz, S. R.; Bur, J. A Three-Dimensional Photonic Crystals Operating at Infrared Wavelengths Nature 1998, 394, 251. [16] Berger, V.; GauthierLafaye, O.; Costard, E. Photonic Band Gaps and Holography J. Appl. Phys. 1997, 82, 60. [17] Scrimgeour, J.; Sharp, D. N.; Blandford, C. F.; Roche, O. M.; Denning, R. G.; Tuberfield, A. J. Three-Dimensional Optical Lithography for Photonic Microstructures Adv. Mater. 2006, 12, 1557. [18] A. Hynninen, H. J. Thijssen, C. M. Vermolen, M. Dijkstra and A. Blaaderen Selfassembly route for photonic crystals with a bandgap in the visible region, Nature Materals 6, pp. 202- 205, (2007) [19] O. Graydon, Nature’s nanosructures colour wings and stones, Opto Lser Europe, 51, pp.31-36 June 1998 [20] R. C. McPhedran, N. A. Nicorovici, D. R. McKenzie, L. C. Botten, A. R. Parker and G. W. Rouse, The Sea Mouse and the Photonic Crystal, Aust. J. Chem. 54, 241-244 (2001). [21] Sanders, J.V.; Sanders, J. V.; Segnit, E. R. (1964). "Structure of Opal". Nature 204: 1151 CHAPITRE II Etude théorique des cristaux photoniques CHAPITRE II : Etude théorique des cristaux photoniques 32 II.1 Introduction L'étude des cristaux photoniques et leurs propriétés spécifiques, mène naturellement à l'étude du comportement de la lumière dans les matériaux à bande interdite photonique. Ces structures périodiques sont régies par les équations de Maxwell. C’est un ensemble de quatre équations différentielles vectorielles qui permettent de modéliser les relations entre les charges, leurs déplacements et les champs électriques et magnétiques. Grâce à l’analogie formelle qui existe entre les équations de Maxwell régissant la propagation des ondes électromagnétiques dans un milieu diélectrique et l’équation de Schrödinger pour les électrons [1], on peut traiter les cristaux photoniques avec les outils et les concepts développés en physique du solide en employant les méthodes de la mécanique quantique et le théorème de Bloch [2]. Cette analogie vient de la périodicité géométrique du cristal. En effet, la périodicité de la constante diélectrique dans l’équation de Maxwell est analogue à la périodicité du potentiel atomique cristallin. Cela nous permet de penser qu’une variation périodique de la permittivité peut conduire à l’apparition de bandes d’énergie interdites pour les photons. II.2 Equations macroscopiques de Maxwell De façon générale, la propagation des ondes électromagnétiques dans un milieu de constante diélectrique ε (r) , y compris la propagation de la lumière dans un cristal photonique, est décrite par les quatre équations de Maxwell (dans le système S.I.): Equation de Maxwell-Faraday 0 B E t ∂ ∇× + = ∂ (2.1) Equation de Maxwell-Ampère D H J t ∂ ∇× − = ∂ (2.2) CHAPITRE II : Etude théorique des cristaux photoniques 33 Equation de Maxwell-Gauss 0 B ∇⋅ = (2.3) Equation de conservation du flux magnétique D ρ ∇⋅ = (2.4) où E désigne le champ électrique, B la densité du flux magnétique, D la densité du déplacement électrique ou la densité du flux électrique, H le champ magnétique, J la densité de courant, ρ la densité de charge électrique, et ∇ l’operateur « nabla » : x y z x y z ∂ ∂ ∂ ∇ = + + ∂ ∂ ∂ (dans les coordonnées cartésiennes). Dans la situation d'un milieu mixte composé de régions de matériau diélectrique homogène qui ne comporte ni charges libres, ni courants libres, dans laquelle la structure ne varie pas avec le temps, nous pouvons mettre 0 ρ = et 0 J = . Généralement, les composantes i D de la densité du déplacement électrique D sont liées aux composantes i E du champ électrique E par une série de puissance [3]: 3 0 , / ( ) i ij j ijk j k j j k D E E E O E ε ε χ = + + ∑ ∑ (2.5) Cependant, pour de nombreux matériaux diélectriques, il est raisonnable d'utiliser les approximations suivantes : • Les champs sont assez faibles pour pouvoir négliger les termes ijk χ (et tous les termes d'ordre supérieur) dans la série (2.5) et observer une relation linéaire entre D et E . • Les matériaux sont macroscopiques et isotropes, de sorte qu’on puisse utiliser une grandeur scalaire pour la constante diélectrique. • La constante diélectrique est supposée indépendante de la fréquence, du moins dans la gamme de fréquences qui nous intéresse pour le système considéré. • On s’intéresse uniquement à des matériaux diélectriques à faibles pertes, ce qui signifie que la constante diélectrique est purement réelle. CHAPITRE II : Etude théorique des cristaux photoniques 34 • Enfin, on suppose la perméabilité magnétique ( ) r µ égale à 1 (ce qui est très proche de la réalité pour la plupart des matériaux diélectriques auxquels on s’intéresse généralement). En considérant les approximations précédentes, on obtient les relations suivantes : 0 ( ) ( ) ( ) D r r E r ε ε = (2.6) 0 ( ) ( ) B r H r ε = (2.7) Avec toutes ces hypothèses, les équations de Maxwell (2.1) – (2.4) deviennent : 0 0 ( , ) ( , ) H r t E r t t µ ∂ ∇× + = ∂ (2.8) 0 ( , ) ( , ) ( ) 0 E r t H r t r t ε ε ∂ ∇× − = ∂ (2.9) 0 ( , ) H r t ∇⋅ = (2.10) 0 ( ) ( , ) r E r t ε ( ∇⋅ = ¸ ¸ (2.11) Comme les équations de Maxwell sont linéaires, il est possible de séparer la dépendance temporelle de la dépendance spatiale et chercher des solutions de type harmonique telles que : ( , ) ( ) i t H r t H r e ω − = (2.12) ( , ) ( ) i t E r t E r e ω − = (2.13) L’insertion des modes harmoniques ci-dessus dans les équations de Maxwell (2.8) - (2.11) donne les deux relations : 0 0 ( ) ( ) E r i H r ωµ ∇× − = (2.14) 0 0 ( ) ( ) ( ) H r i r E r ωε ε ∇× + = (2.15) et conduit aux conditions suivantes : CHAPITRE II : Etude théorique des cristaux photoniques 35 0 ( ) H r ∇⋅ = (2.16) 0 ( ) ( ) r E r ε ( ∇⋅ = ¸ ¸ (2.17) qui ont une interprétation physique simple: il n'ya pas de sources ponctuelles ou des puits des champs de déplacement et magnétique dans le milieu, et les champs électromagnétiques sont des ondes transversales. En partant des équations (2.14) et (2.15), et en employant la relation 0 0 1/ c ε µ = , on peut éliminer l’un des deux champs et obtenir une équation aux valeurs propres pour l’autre : 2 1 ( ) ( ) ( ) H r H r r c ω ε | | | | ∇× ∇× = | | \ . \ . (2.18) Avec les deux équations de divergence (2.16) et (2.17), cette équation nous fournit toutes les informations sur le comportement de ( ) H r . Expérimentalement, on cherche à déterminer le champ électrique au lieu du champ magnétique, car le premier est facilement mesurable. Mais théoriquement la résolution de l’équation transverse électrique devient une tâche plus délicate (la propriété d’hermiticité manque dans cette équation) [4]. Pour cela, il est commode d’éliminer ( ) E r et de conserver l’équation aux valeurs propres pour le champ magnétique ( ) H r . Ensuite, on utilise l’équation (2.15) pour récupérer ( ) E r : 0 ( ) ( ) ( ) i E r H r r ωε ε = (2.19) II.3 Analogie Schrödinger-Maxwell Un photon qui se propage dans un cristal photonique est l'équivalent d'un électron dans un semi-conducteur [5, 6, 7]. Cette analogie électron photon découle de la similitude entre l'équation de Schrödinger régissant la propagation des électrons dans un matériau caractérisé par un potentiel électrostatique périodique et les équations de Maxwell utilisées pour décrire la propagation d'une onde électromagnétique dans un matériau caractérisé par sa constante diélectrique périodique. CHAPITRE II : Etude théorique des cristaux photoniques 36 L'équation de Schrödinger en régime stationnaire pour la fonction d'onde d'un électron dans un potentiel V s'écrit [8] : ( ) 2 2 2 ( ) ( ) ( ) m r U V r r ∇ Ψ = − − Ψ (2.20) ou U est l'énergie de l'électron, m sa masse. Nous avons vu qu'en régime linéaire, l'équation de propagation d'une onde électromagnétique monochromatique dans un matériau avec ( ) r ε était : 2 ( ) ( ) ( ) H r r H r c ω ε | | ( ∇× ∇× = | ¸ ¸ \ . (2.21) Dans ce cas, l'équation de la fonction d'onde d'un électron de masse m dans un potentiel V (équation 2.20) est analogue à l'équation d'onde électromagnétique dans un milieu diélectrique ( ) r ε (équation 2.21). Les équations (2.20) et (2.21) sont deux équations aux valeurs propres. L'équation (2.21) définit les valeurs possibles de la fréquence d'une onde se propageant dans le matériau en l'absence d'excitation extérieure et les amplitudes des champs associés. L’équation (2.20) définit les valeurs possibles de l'énergie d'un électron se propageant librement dans un potentiel et les fonctions d'onde associées. L'énergie de l'électron et la fréquence de l'onde électromagnétique sont les valeurs propres, dictées respectivement par le potentiel et la constante diélectrique. De cette similitude, découlent des propriétés analogues pour les deux systèmes. En identifiant le membre gauche de l'équation maîtresse (2.18) comme un opérateur ˆ Θagissant sur ( ) H r , on arrive à : 2 1 ˆ ˆ ( ) ( ) , ( ) H r H r c r ω ε | | | | Θ = Θ = ∇× ∇× | | \ . \ . (2.22) Dans cette équation, l’opérateur ˆ Θest linéaire et Hermitien (Opérateur dont les éléments de matrice symétriques sont conjugués sur un espace vectoriel complexe). En faisant CHAPITRE II : Etude théorique des cristaux photoniques 37 une intégration par partie deux fois de ˆ ( , ) F G Θ , on peut montrer que ˆ Θest Hermitien pour tous champs de vecteurs ( ) F r et ( ) G r : ( ) 3 * * 3 * 3 1 ˆ ( , ) ( ) 1 1 ˆ ( , ) F G d rF G d r F G d r F G F G ε ε ε Θ = ⋅ ∇× ∇× = ∇× ⋅ ∇× ( | | = ∇× ∇× ⋅ = Θ | ( \ . ¸ ¸ ∫ ∫ ∫ (2.23) En effectuant les intégrations par parties, les termes de surface qui impliquent les valeurs des champs à la limite de l'intégration ont été négligées (les champs qui sont périodiques dans la région de l'intégration, ou les champs qui tendent vers zéro à grandes distances toujours font disparaître les termes de surface). II.3.1 Propriétés des modes harmoniques Les fonctions propres et les valeurs propres de ˆ Θ possèdent des propriétés importantes. Notamment, on va montrer que les valeurs propres sont toujours réelles, et que les fonctions propres associées sont orthogonales [9]. Supposons que ( ) H r est un vecteur propre de ˆ Θavec la valeur propre ( ) 2 / c ω . Le produit intérieur de l’équation maitresse (2.18) avec ( ) H r donne : ( ) ( ) ( ) 2 2 2 2 * * 2 2 ˆ ( ) / ( ) ˆ ( , ) / ( , ) ˆ ( , ) / ( , ) H r c H r H H c H H H H c H H ω ω ω Θ = ⇒ Θ = ⇒ Θ = (2.24) L’opérateur ˆ Θ est Hermitique, donc on peut écrire: ˆ ˆ ( , ) ( , ) H H H H Θ = Θ . En outre, de la définition du produit intérieur ; nous avons la relation : * ˆ ˆ ( , ) ( , ) H H H H Ξ = Ξ pour tout opérateur ˆ Ξ . En exploitant ces deux informations, on obtient : * 2 2 * 2 2 2 2 * 2 2 2 * 2 ˆ ˆ ( , ) ( / ) ( , ) ( , ) ( / )( , ) ( / ) ( / ) ( ) H H c H H H H c H H c c ω ω ω ω ω ω Θ = = Θ = ⇒ = ⇒ = (2.25) CHAPITRE II : Etude théorique des cristaux photoniques 38 et de (2.24) : 2 2 3 1 ˆ ( , ) ( , ) H H H H d r H c ω ε | | = Θ = ∇× | \ . ∫ (2.26) avec ( ) 0 r ε > Par conséquent, toutes les valeurs propres 2 ω sont non-négatives, et ω est réelle. De plus, l’hermiticité de ˆ Θ exige que pour une paire de modes harmoniques quelconques 1 ( ) H r et 2 ( ) H r avec des fréquences différentes 1 ω et 2 ω , le produit intérieur 1 2 ( , ) H H soit nul. Considérons deux modes normalisés 1 ( ) H r et 2 ( ) H r avec des fréquences 1 ω et 2 ω : 2 2 2 2 1 2 1 2 1 2 1 2 2 1 2 2 1 2 2 1 ˆ ˆ ( , ) ( , ) ( , ) ( , ) ( )( , ) 0 H H c H H c H H H H H H ω ω ω ω = Θ = Θ = ⇒ − = (2.27) Si 1 2 ω ω ≠ , alors nous devons avoir 1 2 ( , ) 0 H H = et nous disons que 1 H et 2 H sont des modes orthogonaux. Les modes dégénérés ( 1 2 ω ω = ) ne sont pas nécessairement orthogonaux, mais en tenant compte de la linéarité de l’opérateur ˆ Θ, on peut utiliser des combinaisons linéaires des modes dégénérés qui sont orthogonales. II.3.2 Loi d'échelle Une propriété importante de l'électromagnétisme dans les systèmes diélectriques macroscopiques est qu'il n'existe pas de longueur fondamentale comme pour le rayon de Bohr en physique atomique [10]. On montre en effet que l'équation de propagation (2.18) devient, lorsque la structure étudiée a ses dimensions divisées par un facteur s : 2 1 / ( / ) ( / ) ( / ) s H r s H r s r s c ω ε | | | | ∇× ∇× = | | \ . \ . (2.28) Ainsi, les propriétés à la fréquence ω d'un cristal de constante diélectrique ( ) r ε dont on néglige la dispersion spectrale, sont les mêmes que celles d'un cristal de constante diélectrique ( / ) r s ε à la fréquence / s ω . La géométrie d'un arrangement de matériaux destiné à un travail CHAPITRE II : Etude théorique des cristaux photoniques 39 dans le domaine des micro-ondes peut donc être reprise pour un travail dans le domaine visible : les structures millimétriques servent alors de démonstrateurs avant la réalisation souvent plus lourde d'objets de taille submicronique [11]. II.3.3 Différences et similarités L’équation de propagation électromagnétique est vectorielle alors que celle de Schrödinger est scalaire. Il s’agit d’équations linéaires aux dérivées partielles du deuxième ordre. En ce qui concerne la dérivée temporelle de l’équation de Schrödinger, elle est limitée à l’ordre 1 alors qu’elle atteint l’ordre 2 pour l’équation de propagation de Maxwell. Les électrons sont des fermions. Ils ont un spin demi-entier et suivent la loi de répartition de Fermi. Les photons sont des bosons, ils suivent la loi de répartition de Bose-Einstein. Dans ce cas, il peut y avoir plusieurs particules dans le même état quantique et les bosons tendent naturellement à se regrouper dans le même état. Leur spin est entier. Les photons n’interagissent pas entre eux et leur énergie ne peut pas être modifiée. Ils peuvent être absorbés ou émis, sinon ils conservent leurs fréquences. Electron (Schrödinger) Photon (Maxwell) Périodicité Puits de potentiel électrique carré périodique Constant diélectrique périodique Champ ( , ) ( ) exp( / ) r t r iEt ψ ψ = − ( , ) ( ) exp( ) H r t H r i t ω = Grandeur caractéristique ( ) V r ( ) E r Opérateur Hermitien 2 2 ( ) 2 H V r m − ∇ = + 1 ˆ ( ) r ε | | Θ = ∇× ∇× | \ . Equation aux valeurs propres H E ψ ψ = 2 ˆ ( ) ( ) H r H r c ω | | Θ = | \ . Tableau 2.1 : Récapitulatif Analogie Maxwell-Schrödinger. CHAPITRE II : Etude théorique des cristaux photoniques 40 II.4 Théorème de Bloch Le théorème de Bloch [12], stipule que dans un potentiel périodique, toutes les solutions de l’équation de Schrödinger sont des fonctions dites de Bloch, c’est-à-dire qu’il existe un vecteur k permettant d’écrire : ( ) ( ) ik r k k r u r e ⋅ Ψ = où ( ) k u r est une fonction périodique avec les mêmes périodes que le potentiel. Les fonctions d’ondes des électrons dans un cristal parfait (périodique, infini, sans défaut...) sont donc simplement le produit entre une onde plane et une fonction périodique. L’intérêt de ce théorème est qu'il montre que l'on a uniquement besoin de connaître Ψ sur la maille élémentaire du cristal, les valeurs se reproduisant dans les autres mailles. Les vecteurs k sont appelés vecteurs de Bloch et les fonctions d'ondes sont les fonctions de Bloch. II.4.1 Réseau direct et réseau réciproque Le réseau cristallin (réseau direct) est déterminé par la cellule unitaire. La plus petite cellule unitaire est appelée primitive. elle est sous-tendue par les trois vecteurs fondamentaux 1 a , 2 a , 3 a de telle façon que chaque vecteur de transition du réseau peut être mis sous la forme d'une combinaison linéaire des vecteurs de base : 1 1 2 2 3 3 R n a n a n a = + + (2.29) avec n, n, n : entiers Le réseau réciproque est un réseau de l’espace de Fourier lié au cristal dans lequel le vecteur G , appelé vecteur du réseau réciproque, est un vecteur de translation par lequel l’ensemble du réseau réciproque est construit. G est défini par : 1 1 2 2 3 3 G u b u b u b = + + (2.30) où 1 u , 2 u et 3 u sont des entiers arbitraires et : 2 3 1 1 2 3 2 a a b a a a π × = ⋅ × , 3 1 2 1 2 3 2 a a b a a a π × = ⋅ × , 1 2 3 1 2 3 2 a a b a a a π × = ⋅ × (2.31) CHAPITRE II : Etude théorique des cristaux photoniques 41 Le réseau réciproque, et notamment la première zone de Brillouin, a une grande importance dans la propagation d’onde car les vecteurs d’onde sont toujours tracés dans l’espace de Fourier. II.4.2 Zones de Brillouin Les zones de Brillouin sont des régions qui partitionnent l’espace réciproque associé au cristal. Il en existe une infinité. Elles peuvent être définies à l’aide des plans de Bragg qui sont les plans médiateurs de l’ensemble des vecteurs formés par des combinaisons linéaires des vecteurs du réseau réciproque. La définition de la n-ième zone de Brillouin est la suivante : ensemble des points pouvant être atteint depuis l'origine en croisant n −1 plans de Bragg (Fig. 2.1). Figure 2.1 : Premières zones de Brillouin (ZB) d'un réseau carré. Les plans de Bragg sont tracés avec différentes couleurs. II.4.3 Zone de Brillouin irréductible On appelle « zone de Brillouin irréductible » la plus petite surface qui permet de déduire la relation de dispersion dans tout l’espace réciproque. Elle correspond à la plus petite surface qui peut être utilisée pour reconstruire la première zone de Brillouin (ZB) en utilisant les symétries du réseau réciproque. Pour construire cette zone, nous plaçons au centre de la cellule d’origine G du réseau réciproque pour tracer des vecteurs joignant l’origine aux nœuds voisins de ce même réseau. CHAPITRE II : Etude théorique des cristaux photoniques 42 Nous construisons ensuite les médiatrices de ces vecteurs. La plus petite aire interceptée par ces médiatrices est la zone de Brillouin irréductible. La figure 1-8 donne une représentation graphique des réseaux réciproques, de la première ZB et de la ZB irréductible pour les réseaux 2D carrés et hexagonal. Figure 2.2 : Réseau réel, réciproque, première zone de Brillouin et zone de Brillouin irréductible pour les réseaux 2D (a) carré et (b) triangulaire. II.5 Diagramme de bandes L’analogie avec la physique du solide permet de réutiliser tous les outils de la cristallographie liés à la périodicité du réseau. On peut associer à un cristal photonique une fonction diélectrique périodique ( ) ( ) r r R ε ε = + . Suivant la dimensionnalité de la structure, la constante diélectrique est une fonction périodique du système suivant N =1, 2,3 directions de l’espace, et est invariante selon les (3- N) autres directions. Le vecteur R est une combinaison linéaire des vecteurs de base du réseau direct { i a }. 1 2 3 R la ma na = + + avec l, m, n : entiers (2.32) Dans ce cas, le théorème de Bloch pour un problème aux valeurs propres nous permet de mettre les solutions de l’équation (2.18) sous la forme: ( ) ( ) ik r k k H r u r e ⋅ = (2.33) CHAPITRE II : Etude théorique des cristaux photoniques 43 où ( ) k u r est une fonction périodique, de même période que la structure, et qui est complètement définie par les valeurs qu'elle prend dans la cellule unité : ( ) ( ) k k u r u r R = + . En remplaçant ( ) k H r par sa forme d'onde de Bloch (2.33) dans l'équation maîtresse (2.18), on obtient : 2 2 2 2 ( ) ˆ 1 ( ) ( ) ( ) ( ) 1 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ˆ ( ) ( ) k k ik r ik r k k k k k k k k H H c k e u r e u r r c k ik ik u r u r r c k u r u r c ω ω ε ω ε ω ⋅ ⋅ | | Θ = | \ . | | ∇× ∇× = | \ . | | +∇ × +∇ × = | \ . | | Θ = | \ . (2.34) ˆ k Θ est le nouvel opérateur hermitien qui dépend du vecteur d’onde k : 1 ˆ ( ) ( ) ( ) k ik ik r ε Θ = +∇ × +∇ × (2.35) La résolution d’une telle équation, pour un vecteur d’onde k donné, conduit à un ensemble discret de valeurs propres ( ) n k λ , fonctions du vecteur k , discriminées par un indice de bande entier n. Ces valeurs propres sont reliées aux fréquences propres du cristal par : 2 2 ( ) ( ) n n k k c ω λ = (2.36) C’est l’ensemble des courbes de dispersion des fréquences propres ( ) n k ω qui constitue la structure de bandes du cristal photonique étudié. Ce diagramme de bandes est un élément crucial, car il donne une « cartographie » de tous les états électromagnétiques possibles pouvant exister dans la structure photonique. Les états propres associés à des valeurs propres différentes sont orthogonaux entre eux. A chaque état propre, , n k H correspond une distribution précise du champ électromagnétique obéissant à certaines règles de symétrie. CHAPITRE II : Etude théorique des cristaux photoniques 44 Figure 2.3: Structure de bande d’un réseau carré bidimensionnel de tiges cylindriques diélectriques (ε =8.9); modes TM et TE. Figure 24 : Structure de bande d’un réseau diamant de sphères d’air dans un diélectrique à haute permittivité (ε =13). II.6 Carte des bandes interdites Le calcul du diagramme de bande vu précédemment nous renseigne, entre autres, sur les propriétés (la position et la largeur) des bandes interdites photoniques pour chaque polarisation. Néanmoins, deux paramètres peuvent encore être ajustés afin de jouer sur celles ci : l’indice de la matrice et le facteur de remplissage en air (ou le rapport r/a). Pour un matériau donné (donc pour un indice de la matrice donné), il est intéressant de connaître l’influence du facteur r/a sur la position et la largeur des gaps photoniques : c’est la carte des bandes interdites. CHAPITRE II : Etude théorique des cristaux photoniques 45 La figure 2.5 représente les différentes bandes interdites en fréquence normalisée ωa/2πc (ou a/λ) en modes TE et TM en fonction du rapport r/a dans un cristal photonique de niobate de lithium en configuration carrée. Figure 2.5 : Carte des bandes interdites d’une structure carrée en mode TE et TM pour le niobate de lithium. II.7 Bandes permises et interdites II.7.1 Etude quantique La périodicité du potentiel électrique est modulée par la répartition régulière des ions positifs. Sa périodicité correspond au pas « d » du réseau cristallin. Cette représentation est donnée par la figure 2.6 [13]. Figure 2.6 : Potentiel électrique dans un cristal unidimensionnel CHAPITRE II : Etude théorique des cristaux photoniques 46 Il faut résoudre l’équation de Schrödinger pour déterminer la fonction d’onde de l’électron en utilisant cette forme de potentiel électrique. L’équation d’onde de Schrödinger s’écrit sous la forme suivante : 2 2 2 2 ( ( )) 0 m E V x x ∂ Ψ + − Ψ = ∂ (2.37) Pour simplifier la résolution, on utilise le modèle de Kronig-Penney représenté par un puits de potentiel carré de hauteur de barrière V 0 Figure 2.7 : Puits de potentiel électrique carré périodique La position des atomes est au centre de chaque puits de potentiel, et pour quitter l’atome, l’électron doit lutter contre la force d’attraction représentée par la barrière de potentiel. L’énergie totale E de la particule est supposée telle que et de largeur de barrière b : 0 0 E V ≤ ≤ . Il faut alors résoudre l’équation d’onde de Schrödinger dans les deux régions suivantes, A < x < B et B < x < C , puis appliquer les conditions de continuité et de périodicité aux interfaces. Sur le chemin A-B, on a 0 < x < a et V = 0 d’où l’équation (2.37) s’écrit : 2 1 1 2 2 ( ) 2 ( ) 0 x m E x x ∂ Ψ + Ψ = ∂ (2.38) Sur le chemin B-C, on a a < x < a + b et V =V 0 2 2 0 2 2 2 ( ) 2 ( ) ( ) 0 x m E V x x ∂ Ψ + − Ψ = ∂ d’où l’équation (2.37) s’écrit : (2.39) CHAPITRE II : Etude théorique des cristaux photoniques 47 (avec E <V 0 1 2 2 ( ) sin( ) cos( ) , mE x A x B x α α α Ψ = + = d’après les hypothèses) Ainsi, la solution de l’équation (2.38) est donnée par : (2.40) et la solution de (2.38) est donnée par : 2 0 2 2 ( ) , = ( ) x x m x Ce De V E β β β − + Ψ = + − (2.41) De plus, Ψ(x) est la fonction d’onde solution relative à un motif du puit s de potentiel périodique. La solution générale de l’équation d’onde est une fonction de Bloch qui s’écrit : ( ) ( ) jkx u x x e − = Ψ (2.42) où k est le vecteur d’onde de la particule qui traduit le déplacement de la particule dans le puits de potentiel périodique. Les conditions de continuité doivent être satisfaites au point B sur les fonctions 1 ( ) x Ψ , 2 ( ) x Ψ est sur leurs dérivées ' 1 ( ) x Ψ , ' 2 ( ) x Ψ . Nous obtenons ainsi : 1 2 ( ) ( ) x a x a Ψ = = Ψ = (2.43) ' ' 1 2 ( ) ( ) x a x a Ψ = = Ψ = (2.44) De plus, les conditions de périodicité sont appliquées. Ainsi, la fonction d’onde u(x) (relation 2.34) doit être identique en x = A et en x = C et impose d’avoir : ( 0 )( ) u x u x a b = = = + (2.45) ( 0) ( ) u x u x a b ′ ′ = = = + (2.46) CHAPITRE II : Etude théorique des cristaux photoniques 48 Il est donc nécessaire de résoudre le système suivant : (2.47) afin de déterminer les coefficients A, B, C et D. La résolution de ce système aboutit à la condition suivante : (2.48) Figure 2.8 : Diagramme de dispersion Sur la figure 2.8, apparaît en pointillés la parabole d’équation 2 2 0 ( ) 2 E k k m = . Cette parabole est obtenue dans l’équation (2.48) lorsque b tend vers zéro, c’est-à-dire quand la barrière de potentielle est transparente. Alors, on considère que les électrons sont dans un volume infini. La courbe E(k) ne s’éloigne de cette parabole qu’au voisinage des valeurs k n a π = , en créant chaque fois une bande d’énergie interdite « band gap ». Cette derniere a 1 2 1 2 1 2 ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) sin( ) cos( ) ( ) ( ) ( ) x x jkx x a x a x a x a x x u x o u x a b u u x o x a b x x x A x B x x Ce De u x x e β β α α − + − Ψ = = Ψ = ¦ ¦ ∂Ψ ∂Ψ ¦ = = = ∂ ∂ ¦ ¦ = = = + ¦ ∂ ∂ ¦ = = = + ´ ∂ ∂ ¦ Ψ = + ¦ ¦ Ψ = + ¦ ¦ = Ψ ¦ ¹ | | 2 2 cos( ) ( ) sin( ) ( ) cos ( ) 2 a ch b a sh b k a b β α α β α β αβ − + = + CHAPITRE II : Etude théorique des cristaux photoniques 49 la largeur E G k a π = ± pour . Le domaine de l’espace des k compris entre a π − et a π + s’appelle la première zone de Brillouin de la chaîne atomique. Le domaine constitué des deux segments 2 ; a a π π ( − − ( ¸ ¸ et 2 ; a a π π ( + + ( ¸ ¸ constitue la seconde zone de Brillouin, etc. En utilisant la simplification de Kronig et Penney, l’analyse de l’équation (2.48) est facilité. Pour cela, on suppose que la largeur b de barrière de potentiel tend vers zéro tout en gardant le produit b.V 0 ( ) ch b β constant. Ainsi, tend vers 1, ( ) sh b b β β tend vers 1, 2 b α tend vers 0 et 2 b β tend vers 0 2 2 . mV Dans ce cas précis, l’allure du résultat est inchangée et l’équation ( 2.48) est ainsi simplifiée sous la forme : 0 2 sin( ) cos( ) cos( ) mV a a a ka a α α α + = (2.49) En posant X a α = et 0 2 mV a P = , on obtient alors une réponse du type : sin cos cos( ) X X P ka X + = avec X >0 et P >0 (2.50) En faisant une simple représentation graphique, en prenant une valeur arbitraire pour P (P=4), il est possible de montrer que le membre de gauche de la relation (2.50) admet des valeurs de cette fonction supérieures à 1 et inférieures à -1. Or l’équation (2.50) n’admettra des solutions que lorsque le membre de gauche sera compris entre +1 et -1. Ce comportement met en évidence la notion de « bandes interdites » et de « bandes permises » décrites par la figure 2.9. CHAPITRE II : Etude théorique des cristaux photoniques 50 Figure 2.9 : Bandes permises et bandes interdites II.7.2 Etude électromagnétique Nous allons à présent montrer la similitude existant entre le calcul des modes de propagation électromagnétique dans un matériau périodique et la résolution de l’équation de Schrödinger pour une particule dans un puits de potentiel périodique. Mais avant d’étudier les ces similitudes, établissons la relation de Helmholtz scalaire dérivée des équations de Maxwell. En partant des équations (2.14) et (2.15) et en employant la relation 0 0 1/ c ε µ = , on peut éliminer le champ H et obtenir l’équation aux valeurs propres pour E : 2 ( ) ( ) ( ) E r r E r c ω ε | | ∇×∇× = | \ . (2.51) de plus, on a : 2 2 ( ) ( ( )) ( ) ( ) E r E r E r E r ∇×∇× = ∇ ∇⋅ −∇ = −∇ CHAPITRE II : Etude théorique des cristaux photoniques 51 Alors la relation de propagation d’une onde électromagnétique dans un milieu diélectrique s’écrit, après simplification, de la manière suivante : 2 2 2 ( ) ( ) ( ) 0 E r r E r c ω ε ∇ + = (2.52) Dans un système unidimensionnel, l’équation précédente devient : 2 2 2 2 ( ) ( ) ( ) 0 E x x E x x c ω ε ∂ + = ∂ (2.53) Par analogie entre l’équation différentielle (2.53) et celle de Schrödinger (2.37), il est possible d’identifier le champ électrique E à la fonction d’onde Ψ, et le terme 2 2 ( ) x c ω ε au terme 2 2 ( ( )) m E V x − . En prenant comme hypothèse que la permittivité ε (x) est périodique de période d, on peut mettre de nouveau en évidence la notion de bandes permises et de bandes interdites. Figure 2.10 : Constant diélectrique périodique En reprenant la démarche vue au paragraphe précédent pour le calcul de bandes pour un électron dans un puits de potentiel périodique, nous pouvons faire la résolution de l’équation (2.53) dans les régions A-B et B-C (Figure 2.10). CHAPITRE II : Etude théorique des cristaux photoniques 52 Si 0<x<a alors ( ) 1 a x ε ε = = et : 2 2 2 2 ( ) ( ) 0 AB AB E x E x x c ω ∂ + = ∂ (2.54) Si a < x < a + b alors ( ) b r x ε ε ε = = et : 2 2 2 2 ( ) ( ) 0 BC r BC E x E x x c ω ε ∂ + = ∂ (2.55) Les solutions des équations différentielles (2.54) et (2.55) sont donc respectivement : ( ) sin( ) cos( ) ( ) sin( ) cos( ) AB BC E x A x B x E x C x D x α α β β = + = + (2.56) avec c ω α = et r c ω β ε = En utilisant les faits, qu’au point B la fonction E(x) et sa dérivée E′(x) sont continues, qu’au point C la fonction ε (x) et sa dérivée ε′(x) périodiques, et que le champ électrique E(x) qui vérifie l’équation (2.53) admet une fonction de Bloch pour solution ( ) ( ) jkx k E x x e µ = (où ( ) k x µ est une fonction périodique de même période d que la distribution de permittivité ( ) ( ) k k x d x µ µ + = ), il est possible de montrer que la solution à l’équation (2.53) est : 1 cos( ) cos( ) sin( ) sin( ) cos[ ( )] 2 r r a b a b k a b ε α β α β ε + − = + (2.57) Cette équation présente comme dans le cas des semi-conducteurs, la particularité de n’avoir de solution que lorsqu’elle est comprise entre -1 et +1. Le membre de gauche de l’égalité (2.57) peut être supérieur à +1 ou inférieur à -1. Dans ce cas, il n’y a pas de vecteur d’onde k qui vérifie la relation de dispersion (2.57), donc aucune onde électromagnétique ne se propagera. Comme α et β dépendent tous deux de la pulsation ω, on parle alors de bandes de fréquences interdites. Donc le matériau périodique unidimensionnel empêche les ondes électromagnétiques de se propager à ces fréquences considérées. CHAPITRE II : Etude théorique des cristaux photoniques 53 Figure 2.11 : Diagramme de dispersion 1d Ce paragraphe met en évidence que la périodicité de la permittivité diélectrique peut interdire la propagation des ondes sur une certaine bande de fréquences. Cette notion de périodicité peut être étendue à 2 ou 3 dimensions, mais la nature vectorielle de l’équation de propagation complique considérablement la résolution théorique du problème. II.8 Vitesse de groupe et vitesse de phase La vitesse de groupe est la vitesse de propagation d’un paquet d’ondes. Elle est souvent confondue avec la vitesse de transmission de l’énergie dans le milieu. Elle a un rôle très important dans la propagation de la lumière et les réponses optiques dans les cristaux photoniques. Par conséquent, la connaissance de la vitesse de groupe est essentielle pour comprendre leurs propriétés optiques ; la vitesse de groupe est déterminée par la pente de la courbe de dispersion. La vitesse de phase est la vitesse de déplacement d’une surface équiphase. Pour une phase s’exprimant par : ( , ) x t t kx ϕ ω φ = − + , la vitesse de phase est donc : dx V dt k φ ω = = (2.58) CHAPITRE II : Etude théorique des cristaux photoniques 54 On peut relier la vitesse de groupe à la vitesse de phase en passant dans l’espace des longueurs d’ondes. On part des équations kV φ ω = et : 0 g k k V k ω = ∂ | | = | ∂ \ . (2.59) et on aboutit à la formule de Rayleigh [14] : g V V V φ φ λ λ ∂ = − ∂ (2.60) La vitesse de groupe est égale à la vitesse de phase quand la fréquence angulaire est proportionnelle au vecteur d’onde, c’est-à-dire lorsqu’il n’y a pas de dispersion. C’est le cas d’un système homogène. À partir de cette formule, nous pouvons constater que la vitesse de groupe peut être positive ou négative et être, en valeur absolue, plus ou moins grande que la vitesse de phase. Elle peut passer très rapidement d'une valeur négative à une valeur positive, en particulier en bord de bande interdite [15]. Cette propriété très singulière a donné naissance à de nombreux phénomènes spectaculaires tels que l'effet superprisme ou l'ultraréfraction. Dans le cas d’un cristal photonique à deux ou trois dimensions, la vitesse de groupe doit être remplacée par la relation vectorielle suivante : 0 ( ) g k k V grad ω = = (2.61) Dans ce cas, la vitesse de groupe est portée par la normale à la surface isofréquence ( ) k ω ω = . Elle n’est pas nécessairement colinéaire à la vitesse de phase. Cela apparaît dans les phénomènes d’anisotropie. CHAPITRE II : Etude théorique des cristaux photoniques 55 II.9 Conclusion Ce chapitre à été consacré à l’étude théorique analytique des cristaux photoniques. À partir des équations de Maxwell pour les milieux diélectriques, nous avons établi une équation aux valeurs propres (pour le champ magnétique) similaire à celle de Schrödinger décrivant le mouvement des électrons dans un potentiel périodique d’un cristal atomique. Cette équation, qui peut fournir toutes les informations sur le comportement de la lumière dans des structures photoniques, est vectorielle alors que l’équation de Schrödinger est scalaire. Nous avons également montré que l’opérateur apparaissant dans cette équation, qui est l’analogue de l’opérateur Hamiltonien dans l’équation de Schrödinger, est hermitien, les valeurs propres sont toujours réelles et les fonctions propres associées sont orthogonaux. Nous avons exposé l’analogie entre les semi-conducteurs électroniques caractérisés par des périodicités atomiques et les photons piégés dans des structures photoniques. En exploitant cette analogie formelle, nous avons utilisé les outils et les concepts développés en physique du solide (notamment le théorème de Bloch) et en mécanique quantique pour résoudre le problème aux valeurs propres. La résolution de ce problème donne les modes propres supportés par les structures photoniques et les relations de dispersion qui constituent les structures de bandes photoniques. CHAPITRE II : Etude théorique des cristaux photoniques 56 Bibliographie [1] G. Malpuech, A. Kavokin, G. Panzarini, and A. Di Carlo, Theory of photon Bloch oscillations in photonic crystals, Physical Review B 63, 035108 (2001). [2] J.D. Joannopolous, R.D. Meade, J.N .Winn, Photonic Crystals, Molding the Flow of light Princeton, NJ Princeton University Press (2008). [3] Bloembergen, N. 1965. Nonlinear Optics. NewYork: W. A. Benjamin. [4] H. S. Sözüer, J. W. Haus et R. Inguva, Photonic band : Convergence problems with the plane wave-method, phys. Rev. B 45, p.13962 (1992). [5] E. Yablonovitch, T. J. Gmitter et K. M. Leung, Photonic band structure: The face- centered-cubic case employing non spherical atoms, Phys. Rev. Lett., vol. 67, 2295 (1991). [6] E. Yablonovitch, T.J. Gmitter, R.D. Meade, A.M. Rappe, K.D. Brommer and J.D. Joannopoulos, Donor and acceptor modes in Photonic band structure, Phys. Rev. Lett., vol. 67, 3380 (1991). [7] E. Özbay, B. Temelkuran, Reflection properties and defect formation in photonic crystals Appl, Phys. Lett., vol. 69, 743(1996). [8] N. W. Ashcroft and N. D. Mermin, Solid state physics, 1a Ed., International Thomson Edition, 1976. [9] S. G. Johnson and J. D. Joannopoulos. Photonic Crystals: The Road from Theory to Practice (Kluwer, Boston, 2002). [10] K. Sakoda, Optical Properties of Photonic Crystals, vol. 80 of Springer series in optical sciences, Springer, Berlin Heidelberg New York, (2001). [11] S. D. Cheng, E. Ozbay, S. Mc Calmont and G. Tuttle, K. –M. Ho, Optimized dipole antennas on photonic band gap crystals, App. Phys. Lett. 67,p. 3399-3401 (1995). [12] C. Kittel, Introduction to Solid State Physics, John Wiley & Sons, Inc. 1996. [13] Yannick MERLE, Etude de la dispersion électronique dans les matériaux périodiques diélectriques bidimensionnels, Thèse de doctorat n° 47-2003, Université de Limoges, Novembre 2003. [14] J. M. Lourtioz, Photonic Crystals: Towards Nanoscale Photonic Devices. New York: Springer, 2008. [15] B. Gralak, S. Enoch and G. Tayeb, Anomalous refractive properties of photonic crystals, Journal of Optical Society of America A 17, 1012-1020 (2000). CHAPITRE III Méthodes de simulation numérique CHAPITRE III : Méthodes de simulation numérique 58 III.1 Introduction La complexité de la fabrication et de la caractérisation des cristaux photoniques aux fréquences optiques rend coûteuses en temps et argent les études expérimentales systématiques sur des dispositifs des cristaux photoniques. Le développement des méthodes de modélisation numérique précises et rapides reste donc primordial pour l’étude de ces structures. Les calculs théoriques des cristaux photoniques sont, en principe, exacts, parce que les équations de Maxwell sont dérivées des premiers principes. Par conséquent, la puissance des calculs est comparable à celle des expériences dans la caractérisation des cristaux photoniques. Parmi les modèles théoriques traitant les cristaux photoniques, on doit distinguer tout d'abord deux catégories qui dépendent de la taille finie ou infinie des structures et puis de leur dimensionnalité (1D, 2D ou 3D). Dans la première catégorie, les méthodes des différences finies ou FDTD, et les méthodes basées sur les matrices de transfert sont le plus souvent utilisées, elles permettent de calculer la réponse spectrale d'un dispositif, ainsi que les propriétés de réflexion et de transmission de la difractions par le réseau du cristal photonique, on peut citer aussi les théories de difractions (RCWA). Les principales techniques utilisées dans la deuxième catégorie traitant des cristaux de taille infinie sont basées sur la décomposition en ondes planes. Avant de détailler les deux méthodes numériques que nous avons choisies pour développer notre logiciel de simulation, nous avons jugé intéressant de faire un tour d’horizon des différents types de méthodes numériques les plus utilisées. CHAPITRE III : Méthodes de simulation numérique 59 III.2 Méthode de décomposition en ondes planes La méthode de décomposition en ondes planes PWE (Plane Wave Expansion) [1] est une méthode de résolution dans le domaine fréquentiel des équations de Maxwell. Elle est basée sur la décomposition en ondes planes du problème électromagnétique. Cette méthode est essentiellement utilisée pour analyser les propriétés dispersives des matériaux à bandes interdites photoniques et permet de déterminer la fréquence, la polarisation, la symétrie et la distribution du champ pour les modes d'une structure photonique [2]. Figure 3.1 : Exemple d’un diagramme de dispersion pour un cristal photonique tridimensionnel. D’un point de vue général, un milieu de constante diélectrique périodique ( ) r ε induit des modifications sur la propagation de l’onde électromagnétique. La recherche des solutions harmoniques ( , ) ( ) j t H r t H r e ω = et ( , ) ( ) j t E r t E r e ω = au système d’équations de Maxwell dans un milieu sans source et non absorbant aboutit à deux équations d’ondes découplées pour le champ électrique (équation 3.1) et magnétique (équation 3.2). ( ) 2 1 ( ) ( ) ( ) E r E r r c ω ε ∇× ∇× = (3.1) 2 1 ( ) ( ) ( ) H r H r r c ω ε ∇× ∇× = (3.2) CHAPITRE III : Méthodes de simulation numérique 60 L’équation d’onde dans un milieu sans perte est formellement analogue à l’équation de Schrödinger décrivant la fonction d’onde d’un électron. La différence est la nature de l’équation d’onde, vectorielle dans le cas des photons, scalaire dans le cas des électrons. Dans les matériaux à bandes interdites photoniques diélectriques, l’information sur la périodicité est contenue dans la fonction diélectrique. On peut écrire pour tout vecteur R du réseau direct la relation suivante : ( ) ( ) r r r R r ε ε + = (3.3) Compte tenu de la périodicité de la permittivité, l’équation d’onde se résout en décomposant ( ) E r et ( ) H r en ondes planes : il s’agit de décomposer le problème en série de Fourier spatiale car le milieu est périodique. La décomposition en série de Fourier impose la définition d’une base de vecteurs orthogonaux aux vecteurs de bases qui portent la périodicité du problème. Les vecteurs de base caractérisent le matériau et forment le réseau direct. A ce dernier, on fait correspondre un réseau réciproque formant la base de décomposition en série de Fourier. Le théorème de Bloch permet d’écrire la permittivité sous la forme suivante : ( ) ( ) ikr k r u r e ε = (3.4) Dans cette dernière relation, ikr e est une onde plane et ( ) k u r est une fonction d’onde appelée souvent « fonction de Bloch ». Les champs électromagnétiques sont décomposés de la même façon. Ceci permet de réduire l’équation faisant apparaître E (ou H ) à une équation aux valeurs propres qui peut être résolue numériquement par les algorithmes éprouvés en physique du solide. Intérêts L'avantage principal de la méthode PWE est sa performance. Pour les structures les plus simples, la vitesse de calcul est très grande et les calculs ne consomment que peu de ressources informatiques. La méthode PWE permet de calculer la structure de bandes non seulement pour les cristaux photoniques les plus simples comme les structures périodiques infinies, mais aussi pour des structures plus compliquées telles que les structures avec défauts ponctuels. CHAPITRE III : Méthodes de simulation numérique 61 En outre, cette méthode permet aussi de calculer la distribution de champ pour chacun des états propres du cristal photonique. Limitations La méthode PWE est incapable de traiter les milieux dissipatifs. L'opérateur différentiel dans l'équation de Helmholtz est hermitien seulement en cas de milieu sans pertes. L’algorithme standard (sans amélioration) de cette méthode ne peut pas traiter les matériaux dispersifs. III.3 Méthode des différences finies dans le domaine temporel La méthode des différences finies dans le domaine temporel « Finite Difference Time Domain, FDTD » a été initialement proposée par Kane S. Yee en 1966 [3] pour résoudre les problèmes impliquant les équations de Maxwell dans les milieux isotropes avec les conditions aux limites. Il a développé les premiers algorithmes de calcul concernant l’approche temporelle. Ensuite, la méthode FDTD à été appliquée aux cristaux photoniques. Sa versatilité permet de simuler la plupart des systèmes [4]. Cette méthode est particulièrement intéressante pour connaître la réponse spectrale d’un système non nécessairement périodique et pour calculer les distributions de champ dans des structures de dimensions finies [5]. La FDTD permet non seulement de calculer les diagrammes de bandes mais aussi de simuler l'évolution temporelle du champ électromagnétique dans les cristaux photoniques, ce qui permet d'avoir des informations sur de nombreuses autres grandeurs comme le vecteur de Poynting ou l'énergie électromagnétique stockée. Un maillage dans l’espace réel est réalisé afin de discrétiser les champs et de pouvoir estimer leurs dérivées. Les champs sont propagés dans le temps de proche en proche à partir d’une certaine distribution de départ donnée par l’utilisateur et pendant un certain temps T. On peut donc observer le régime transitoire du système et voir le régime permanent s’établir. De plus, il est possible de remonter à la réponse fréquentielle de la structure. En effet, connaissant l’évolution temporelle des champs en tous points du système, le calcul du spectre en ces points peut être effectué par transformée de Fourier. Le spectre d’un point du système sera alors divisé par celui de la source. CHAPITRE III : Méthodes de simulation numérique 62 Figure 3.2 : Cellules de Yee à 3 dimensions Les conditions aux limites constituent un point très important dans cette méthode du fait du traitement presque exclusif de structures finies. Il est possible d’utiliser à la fois des conditions périodiques et des conditions absorbantes sur une même structure. Ceci est intéressant, par exemple, pour la simulation d’un guide à cristal photonique. Parmi les conditions d’absorption aux limites, les plus utilisées sont les conditions de Mur [6] ou de PML (Perfectly Matched Layer) [7]. Cette dernière, basée sur le placement d’un absorbant artificiel sur le bord de la cellule, est définie pour avoir une adaptation d’impédance avec le vide et pour éliminer toute réflexion. Le principal inconvénient de la FDTD est qu'elle nécessite de longs temps de calculs et un espace mémoire important, en particulier pour les calculs 3D. C'est pourquoi, nous nous limiterons souvent à des calculs 2D. Intérêts Un code numérique de type FDTD est relativement simple à mettre en œuvre et rapide (le nombre d'opérations arithmétiques mises en œuvre à chaque itération est faible). La méthode est utilisable à la fois dans les domaines temporel et fréquentiel. CHAPITRE III : Méthodes de simulation numérique 63 On peut modéliser des structures aux géométries complexes dont les matériaux peuvent être anisotropes ou inhomogènes. Le schéma FDTD est explicite en temps : il n'y a pas de système linéaire à résoudre et il demande peu de stockage. Ce schéma est naturellement centré, d'où une précision à l'ordre 2 en espace et en temps. Il est aisément parallèlisable : c'est un schéma différences finies, donc local, par opposition aux méthodes intégrales. La méthode FDTD donne la possibilité d’intégrer dans l’algorithme de résolution de nombreuses sources (ondes planes, modes guidés, dipôles oscillants,…). Le calcul de toutes les composantes du champ à tous les instants et sur tous les domaines permet d’obtenir des spectres en fréquence en divers points de la structure grâce à l’utilisation de la transformée de Fourier. Limitations Le schéma de Yee impose une discrétisation en mailles régulières. Cette méthode nécessite l'utilisation d’une grille cartésienne ; cela entraîne des difficultés pour faire des raffinements locaux. En particulier, les géométries courbes doivent être approchées par des courbes en escalier. Néanmoins, un certain nombre de solutions plus ou moins satisfaisantes existent déjà. Il faut des longueurs d’arêtes entre λ /16 et λ /20 pour avoir la convergence en maillage. Il est difficile de prendre en compte des matériaux dispersifs. Les conditions aux limites approchées introduisent une approximation qui peut avoir une influence sur la précision du calcul. FDTD étant un schéma explicite, le pas de temps maximal est relié au pas d'espace par une condition de stabilité. Le pas de maillage et par conséquent le pas de temps sont liés à la plus petite longueur d'onde λ considérée. Un maillage typique aura au moins 10 mailles par longueur d'onde pour garantir une précision acceptable sur la solution. Donc, pour traiter un objet de taille 100λ , il faudra un maillage 1000 x 1000 x 1000, soit 6 milliards d'inconnues scalaires, d'où un coût mémoire très important sur de gros maillages. CHAPITRE III : Méthodes de simulation numérique 64 III.4 Méthode des éléments finis La méthode des éléments finis (Finite Elements Method, FEM) est conceptualisée par A. Hrennikoff et R. Courant dans les années 1940 pour résoudre des problèmes de mécanique de structures [8]. Quelques années plus tard, elle est introduite en électromagnétisme pour résoudre les équations de Maxwell. Cette méthode repose sur un découpage de l'espace selon un maillage. D'habitude, l'on choisit un maillage tétraédrique qui a l’avantage de s’adapter facilement aux structures complexes mais rien n’empêche de le modifier suivant la géométrie du domaine et de mailler plus finement certains endroits que d’autres (Fig. 3.3). La résolution de l’équation aux dérivées partielles sur chaque nœud du réseau des éléments finis donne une solution approximée par des fonctions d’interpolation. Figure 3.3. Exemples de maillages adaptatifs utilisés en FEM. L’équation aux dérivées partielles ne peut pas être résolue directement dans le maillage. Elle doit pour cela être écrite sous une forme variationnelle. Cette forme contient les informations de l’équation aux dérivées partielles et les conditions aux limites. Le principe consiste à minimiser ou maximiser l’énergie autour de la solution réelle. Parallèlement, la solution est approximée sur un élément par des fonctions d’interpolation i α . La même base de fonctions d’interpolation est utilisée pour tous les éléments qui sont très souvent une famille de fonctions polynômes. La somme de ces fonctions pondérées par des coefficients i ψ nous donne la solution sur un élément : 1 1 N élément i i i u ψ α = = ∑ (3.4) CHAPITRE III : Méthodes de simulation numérique 65 Cette méthode est largement employée dans l’étude des cristaux photoniques. Elle permet d’avoir accès aux coefficients de réflexion et de transmission du BIP, à des cartes de champ et aux diagrammes de rayonnements. Enfin, la solution u est remplacée dans la formulation variationnelle par la solution approximée. Un système d’équations dont les variables sont les coefficients i ψ de chaque élément est obtenu et la solution finale en est extraite. Figure 3.4: Les différentes étapes de la méthode FEM Intérêts Les domaines peuvent avoir des tailles très contrastées et des formes quelconques. Les matériaux utilisés peuvent être de nature très différentes : anisotropes, métaux, parfaitement conducteurs, etc. Les matériaux peuvent être dispersifs : leurs paramètres peuvent dépendre de la fréquence de calcul. C’est un avantage propre à toutes les méthodes harmoniques. Limitations Son implémentation est peu aisée et son utilisation demande des ressources de calcul importantes. L’obligation d’utiliser des conditions aux limites absorbantes rajoute non seulement des degrés de liberté au système mais également une approximation supplémentaire. L’équation aux dérivées partielles + Conditions aux limites Formulation variationnelle Fonctions d’interpolation Système d’équations à résoudre Définition de la géométrie Maillage du domaine CHAPITRE III : Méthodes de simulation numérique 66 III.5 Méthode rigoureuse des ondes couplées La méthode rigoureuse des ondes couplées (Rigorous Coupled Wave Analysis, RCWA) est une méthode différentielle de diffraction par les réseaux. Elle a été décrite pour la première fois en 1981 par Moharam et Gaylord pour des réseaux plans diélectriques ou métalliques modulés sinusoïdalement [9]. Elle fut généralisée aux réseaux présentant des reliefs de surface puis au cas tridimensionnel d’incidence conique. La méthode a également été étendue aux cas de réseaux anisotropes. Elle fut rigoureusement formulée un peu plus tard [10]. La méthode RCWA est basée sur la décomposition du champ électromagnétique et de la permittivité diélectrique en séries de Fourier. Pour cette raison, elle est appelée « méthode modale par expansion de Fourier ». Cette méthode est rigoureuse car elle résout les équations de Maxwell dans toute leur généralité sans recourir à des approximations. Intérêts Elle est relativement facile à programmer et à mettre en œuvre sur des géométries simples. Sa convergence est prouvée en fonction du nombre de couches prises en compte dans la modélisation et du nombre de modes de Bloch. Les ressources informatiques nécessaires à un calcul sont en général réduites par rapport à des méthodes d’éléments finis. Limitations Lorsqu’il y a de forts contrastes d’indices entre les matériaux qui constituent la structure modélisée, il faut garder beaucoup de termes dans la décomposition de Fourier de la permittivité (ce qui revient à faire sa transformation de Fourier discrète sur un grand nombre de termes) car il faut prendre en compte correctement les fortes discontinuités aux interfaces entre matériaux. Il est d’une part difficile de savoir quel est le nombre minimum de termes à garder dans la série de Fourier pour avoir une précision acceptable. D’autre part, garder un grand nombre de termes dans la série de Fourier peut être très handicapant en terme de ressources mémoire et de temps de calcul. La convergence du calcul impose de prendre un nombre de couches suffisamment grand pour modéliser correctement la structure réelle. CHAPITRE III : Méthodes de simulation numérique 67 Si la structure étudiée comprend des détails de taille caractéristique très inférieure au pas du réseau, la méthode RCWA est peu précise et très coûteuse. La discrétisation en couches conduit à une discrétisation approchée par des pavés des surfaces. III.6 Méthode de la ligne de transmission La méthode de la ligne de transmission (Transmission Line Matrix, TLM) a été inventée en 1971 par Peter B. Johns associé à Raymond L. Beurle [11], et constitue une approche cousine de l’algorithme FDTD. Elle résout les équations de Maxwell dans le domaine temporel. L'approche de base de la méthode TLM est d’obtenir un modèle discret, résolu ensuite exactement par des moyens numériques; les approximations ne sont introduites qu’au stade de la discrétisation. Cela contraste avec l'approche traditionnelle dans laquelle un modèle continu idéalisé est d'abord obtenu, ensuite résolu par approximations. L’algorithme TLM est basé sur une analogie simple entre la propagation des ondes électromagnétiques dans un milieu et la propagation des tensions et des courants dans un réseau de lignes de transmission. Cette analogie est naturelle si nous considérons la ressemblance frappante des équations de Maxwell et de l’équation du télégraphiste. Il suffit alors de simuler la propagation des tensions et courants dans un réseau de lignes de transmission adéquat pour en déduire celle des ondes électromagnétiques dans un environnement donné. Pour les systèmes électromagnétiques, le modèle discret est formé par le remplissage du domaine de calcul par un réseau de lignes de transmission de manière à ce que la tension et le courant donnent des informations sur les champs électriques et magnétiques. Par conséquent, elle considère le domaine de calcul comme un maillage de ligne de transmissions. Le point où les lignes de transmission se croisent est considéré comme un nœud. A chaque pas de temps, les impulsions de tension arrivent sur chaque nœud de la ligne de transmission. Ces impulsions sont ensuite dispersées pour produire une nouvelle série d'impulsions qui deviennent des nœuds adjacents incidents lors du prochain pas de temps. La relation entre les impulsions incidentes et dispersées est déterminée par la matrice de diffusion, qui doit être compatible avec les équations de Maxwell. La figure 3.5 considère un exemple simple du maillage à deux dimensions de la méthode avec une tension d’impulsion partant du nœud central. Cette impulsion sera partiellement réfléchie et transmise telle que la théorie de la ligne de transmission le décrit. CHAPITRE III : Méthodes de simulation numérique 68 Figure 3.5 : Un exemple de la méthode TLM en 2D: une tension d'impulsion incidente est diffusée deux fois. Le nœud le plus couramment utilisé en 3 dimensions est le nœud symétrique condensé que P.B. Johns créa en 1987 [12]. Il se compose de 12 ports pour que les deux polarisations TE (transverse électrique) et TM (transverse magnétique) des champs soient fixées pour chacune des 6 faces du maillage (Fig. 3.6). D’autres éléments, peuvent être ajoutés au nœud afin que les différentes propriétés du matériau puissent être représentées. Figure 3.6 : Paramètres du nœud condensé symétrique (SCN). CHAPITRE III : Méthodes de simulation numérique 69 III.7 Méthode des matrices de transfert La méthode des matrices de transfert (Transfer Matrix Method, TTM) permet de calculer les coefficients de réflexion et de transmission pour les cristaux photoniques de taille finie. Ces matériaux peuvent être parfaits ou dopés. Elle peut être utilisée pour le calcul de structure de bande pour un cristal parfait. Elle a été adaptée au cas des cristaux photoniques par Pendry au début des années 90 [13] et plus tard par Reynolds [14]. Figure 3.7 : Diagrammes de transmission calculés avec la méthode des matrices de transfert (trait plein) et la méthode des éléments finis. En fait, la méthode TTM permet d’exprimer le champ électromagnétique sur une couche en fonction du champ sur la couche précédente. L’évolution du champ dans le cristal est alors calculée de couche en couche, ce qui permet d’obtenir les coefficients de transmission et de réflexion. L’avantage de cette méthode par rapport à la méthode FDTD c’est qu’elle occupe peu d’espace mémoire. III.8 Approches hybrides Comme nous l'avons vu précédemment, chaque méthode de modélisation est très spécifique. A l'heure actuelle, aucune méthode ne peut répondre pleinement à l'ensemble des problèmes posés, chacune ayant ses avantages et ses inconvénients. CHAPITRE III : Méthodes de simulation numérique 70 On se projette donc de plus en plus vers des approches multi-méthodes, appelées aussi méthodes hybrides. Mais cela n'est pas sans poser de problèmes pour définir des interfaces entre méthodes, même quand celles-ci semblent assez proches dans leur formulation. Dans les approches hybrides, l'espace d'étude est décomposé en différents volumes où l’on recherche à utiliser la méthode d'analyse la plus appropriée. Sur les surfaces constituant les interfaces, des techniques de raccordement sont alors élaborées pour assurer la bonne articulation des méthodes et la continuité des phénomènes électromagnétiques. Si les approches hybrides sont bien appropriées aux méthodes fréquentielles entre elles, elles ont été moins utilisées dans le domaine temporel. On peut noter toutefois des réalisations dans le couplage de méthode FDTD avec la FVTD [15] qui permet un maillage conformé équivalent à la FEM . D’autres codes ont été développés en couplant la FDTD avec la version temporelle de la MoM (MoMTD) [16]. Dans le domaine temporel, on a privilégié les méthodes diakoptiques permettant de scinder le problème en sous problèmes résolus distinctement. On peut noter, par exemple, la FDTD à résolution multiples (MR-FDTD) ou la FDTD à grilles duales (DG-FDTD) [17]. Certains principes fondamentaux de l’électromagnétisme facilitent l’articulation des méthodes entre elles. C’est le cas d’un des principes d’équivalence des champs. On l’illustre ici en considérant deux régions dont l’une contient toutes les sources. Ces dernières créent aux interfaces (surface de Huyguens), des champs dont la connaissance suffit à résoudre le problème dans l’autre région. Les approches les plus simples de ce principe ne prennent pas en compte les effets de couplage qui peuvent exister entre les différentes régions du problème. Ceux-ci peuvent être considérés à partir d’une démarche itérative. En effet, si l’un des objets présent dans la région secondaire entraîne une modification importante des champs dans la région contenant les sources, ceux-ci doivent être recalculés pour bien tenir compte de leur évolution aux interfaces. Le gain de temps de telles méthodes peut être contrarié si le couplage entre sous domaines est fort ou si les interfaces entre ces mêmes domaines sont mal définies. CHAPITRE III : Méthodes de simulation numérique 71 Bibliographie [1] K. M. Ho, C. T. Chan, et M. Soukoulis, Existence of a photonic gap in periodic structures, Phys. Rev. Lett., vol. 65, p. 3152, 1990. [2] S. G. Johnson, S. Fan, P. R. Villeneuve, J. D. Joannopoulos, et L. A. Kolodziejski, Guided modes in photonic crystal slabs, Phys. Rev. B, vol. 60, pp. 5751-5758, 1999. [3] S. K. Yee, Numerical solution of initial boundary value problems involving Maxwell’s equations in isotropic media. IEEE Trans. Antennas and propagation, Vol. 14, pp 302- 307, (1966). [4] Taflove, A. and S.C. Hagness, Computational electrodynamics: the finite-difference time-domain method. 2000: Artech House, Boston. [5] C. T. Chan, Q. L. Yu et K. M. Ho, Order N spectral method for electromagnetic waves, Physical Review B 51, p. 16635 (1995). [6] G. Mur, Absorbing boundary conditions for the finite-difference approximation of the time-domain electromagnetic field equations, IEEE Trans. Electromagnetic Compatibility 23, p. 377 (1981). [7] J.P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves, Journal of Computational Physics 114, p. 185 (1994). [8] R. L. Courant, Variational Methods for the Solution of Problems of Equilibrium and Vibration, Bulletin of the American Mathematical Society 49: 1-23., 1943. [9] M. G. Moharam et T. K. Gaylord, Rigorous coupled-wave analysis of planar grating diffraction, Journal of the Optical Society of America 71, 811 (1981). [10] L. Li et C. W. Haggans, Convergence of the coupled-wave method for metallic lamellar diffraction gratings, Journal of the Optical Society of America A 10, 1184 (1993). [11] P. B. Johns et R. L. Beurle, Numerical solution of 2-dimensional scattering problems using a transmission-line matrix, Proceedings IEE, vol. 118, p. 1203–1208, sept. 1971. [12] P. B. Johns, A symmetrical condensed node for the TLM method, IEEE Transactions on Microwave Theory and Techniques, vol. 35, p. 370–377, avril 1987. [13] J. B. Pendry, Calculating photonic band structure, J. Phys.: Condens. Matter, Vol. 8, 1085, 1996. [14] A. L. Reynolds, Translight Software Manual, University of Glasgow(2000). CHAPITRE III : Méthodes de simulation numérique 72 [15] K. S. Yee and J. S. Chen, The finite-difference time-domain (FDTD) and the finite- volume time-domain (FVTD) methods in solving Maxwell’s equations, IEEE Trans. Antennas Propagat. 45, 354 (1997). [16] A. R. Bretones, R. Mittra and R. G. Martin, “A Hybrid Technique Combining the Method of Moments in the Time Domain and FDTD”, IEEE Microwave and Guided Wave Letters, vol. 8, n°8, August 1998, pp. 281-283. [17] R. Pascaud, R. Gillard, R. Loison et al., Une amélioration de la MR-FDTD basée sur son hybridation avec la DG-FDTD, JNM 2007, Toulouse, Mai 2007. CHAPITRE IV Méthode des différences finies dans le domaine temporel CHAPITRE IV : Méthode des différences finies dans le domaine temporel 74 IV.1 Introduction La méthode des différences finies dans le domaine temporel, plus connue sous l’appellation FDTD (Finite Differences Time Domain), a été formulée pour la première fois en 1966 par Yee [1] et ce n’est qu’en 1975 qu’une série d’articles préconisant ses futures applications sont apparus [2] [3] [4]. Grâce à ses avantages et à l’outil informatique devenant de plus en plus performant, la FDTD n’a cessé de gagner d’utilisateurs pour des applications de plus en plus variées. Le schéma proposé par Yee permet de surmonter la difficulté due à la dépendance simultanée des champs électrique et magnétique entre eux. En effet, on obtient un schéma numérique explicite permettant le calcul du champ électromagnétique dans tout le volume d’étude en fonction du temps. Les composantes des champs électriques et magnétiques sont décalées d’un demi-pas spatial et calculées de manière alternative et itérative, respectivement à des multiples pairs et impairs du demi-pas temporel. Cette méthode peut simuler le comportement d’une onde électromagnétique dans tout type de milieu (diélectrique, métal, plasma....), tout en tenant compte des formes géométriques les plus complexes des objets pouvant constituer le système. Elle trouve ses principales applications dans les domaines de la conception (antennes et circuits), de la compatibilité électromagnétique, de la diffraction, de la propagation et de la dosimétrie électromagnétique (interaction ondes - matière vivante). La FDTD est une méthode bien adaptée pour la modélisation des cristaux photoniques, principalement parce qu’elle permet un accès aux caractéristiques dynamiques des structures (facteur de qualité de résonateur, transmission, réflexion). Elle permet l’utilisation de sources des profils spatiaux (modes guidés, ondes planes, source ponctuelle) et des profils temporels (harmoniques, impulsions) divers, ainsi que l’utilisation de conditions aux limites adaptées aux fortes diffractions qui apparaissent dans les cristaux photoniques. Nous allons exposer dans ce chapitre, la base théorique et les principaux points clés de la méthode FDTD : la discrétisation des équations de Maxwell aux sens des différences finies centrées, le critère de stabilité et la dispersion numérique due au maillage spatial, les conditions aux limites, le traitement des interfaces entre différentes couches diélectriques , le traitement des matériaux dispersifs… CHAPITRE IV : Méthode des différences finies dans le domaine temporel 75 IV.2 Equations de Maxwell dans l’espace cartésien On reprend les équations de Maxwell-Faraday et de Maxwell-Ampère dans le domaine temporel sous leur forme différentielle [5] : B E M t ∂ = −∇× − ∂ (4.1) D H J t ∂ = ∇× − ∂ (4.2) avec J E σ = , * M H σ = , σ est la conductivité électrique, * σ est la résistivité magnétique équivalente. Dans un matériau linéaire, isotrope et non dispersif, les champs B et H d’une part et D et E d’autre part sont reliés par les relations suivantes : 0 r D E E ε ε ε = = ; 0 r B H H µ µ µ = = (4.3) Où ε et µ représentent la permittivité électrique et la perméabilité magnétique, 0 ε est la permittivité électrique du vide, 0 µ est la perméabilité magnétique du vide, r ε est la permittivité relative, r µ est la perméabilité relative. La projection de ces deux équations (4.1) et (4.2) sur un repère cartésien (x, y, z) donne six équations relatives aux différentes composantes des champs électrique et magnétique : 1 1 1 y x z x y z x y y z x z E H E M t z y H E E M t x z E H E M t y x µ µ µ ∂ ( ∂ ∂ = − − ( ∂ ∂ ∂ ¸ ¸ ∂ ∂ ∂ ( = − − ( ∂ ∂ ∂ ¸ ¸ ∂ ( ∂ ∂ = − − ( ∂ ∂ ∂ ¸ ¸ (4.4) CHAPITRE IV : Méthode des différences finies dans le domaine temporel 76 1 1 1 y x z x y x z y y z x z H E H J t y z E H H J t z x H E H J t x y ε ε ε ∂ ( ∂ ∂ = − − ( ∂ ∂ ∂ ¸ ¸ ∂ ∂ ∂ ( = − − ( ∂ ∂ ∂ ¸ ¸ ∂ ( ∂ ∂ = − − ( ∂ ∂ ∂ ¸ ¸ (4.5) IV.3 Réduction à deux dimensions Dans le cas où les cristaux photoniques sont périodiques suivant deux directions (x et y par exemple) et infini suivant la troisième direction (z), on peut distinguer deux types de propagation : propagation dans le plan (in-plane, k z nul) et hors plan (off-plane, k z La propagation in-plane se fait dans le plan et la variation des champs s’annule suivant la troisième direction. Le système d’équations précédent se simplifie et se subdivise en deux sous-systèmes donnant naissance aux deux polarisations transverse électrique : TE non nul). z et transverse magnétique TM z 1 1 1 x z x y z y y z x z H E M t y H E M t x E H E M t y x µ µ µ ( ∂ ∂ = − − ( ∂ ∂ ¸ ¸ ∂ ∂ ( = − ( ∂ ∂ ¸ ¸ ∂ ( ∂ ∂ = − − ( ∂ ∂ ∂ ¸ ¸ . Pour illustrer ce cas, supposons que le cristal est périodique suivant les directions x et y et est infini suivant z. Les dérivées par rapport à z s’annulent. Les équations (4.4) et (4.5) impliquent : (4.6) 1 1 1 x z x y z y y z x z E H J t y E H J t x H E H J t x y ε ε ε ( ∂ ∂ = − ( ∂ ∂ ¸ ¸ ∂ ∂ ( = − − ( ∂ ∂ ¸ ¸ ∂ ( ∂ ∂ = − − ( ∂ ∂ ∂ ¸ ¸ (4.7) Remarquons que l’évolution de E z , H x et H y est indépendante de celle de E x , E y et H z . De ce fait résulte deux systèmes d’équations indépendants: l’un décrivant la polarisation TE z et l’autre, la polarisation TM z . CHAPITRE IV : Méthode des différences finies dans le domaine temporel 77 IV.3.1 Polarisation TE Les composantes électriques sont transverses, elles sont dans le plan de périodicité du cristal photonique. La polarisation TE z 1 1 1 x z x y z y y z x z H E M t y H E M t x H E H J t x y µ µ ε ( ∂ ∂ = − − ( ∂ ∂ ¸ ¸ ∂ ∂ ( = − ( ∂ ∂ ¸ ¸ ∂ ( ∂ ∂ = − − ( ∂ ∂ ∂ ¸ ¸ est définie par le système suivant : (4.8) IV.3.2 Polarisation TM Dans ce cas, ce sont les composantes magnétiques qui sont transverses. Cette polarisation est régie par le système suivant : 1 1 1 x z x y z y y z x z E H J t y E H J t x E H E M t y x ε ε µ ( ∂ ∂ = − ( ∂ ∂ ¸ ¸ ∂ ∂ ( = − − ( ∂ ∂ ¸ ¸ ∂ ( ∂ ∂ = − − ( ∂ ∂ ∂ ¸ ¸ (4.9) IV.3.3 Propagation off-plane La propagation off-plane est caractérisée par une constante de propagation k z ( , , , ) ( , , ) exp( ) ( , , , ) ( , , ) exp( ) z z E x y z t E x y t ik z H x y z t H x y t ik z = = non nulle suivant z. Dans ce cas, les vecteurs champs électriques et magnétiques peuvent s’écrire comme suit : (4.10) CHAPITRE IV : Méthode des différences finies dans le domaine temporel 78 Par conséquent, les dérivées par rapport à z dans le système d’équations de Maxwell (4.4) et (4.5) sont calculées de manière analytique. Ce système devient : 1 1 1 x z z y x y z z x y y z x z H E ik E M t y H E ik E M t x E H E M t y x µ µ µ ( ∂ ∂ = − − ( ∂ ∂ ¸ ¸ ∂ ∂ ( = − − ( ∂ ∂ ¸ ¸ ∂ ( ∂ ∂ = − − ( ∂ ∂ ∂ ¸ ¸ (4.11) 1 1 1 x z z y x y z z x y y z x z E H ik H J t y E H ik H J t x H E H J t x y ε ε ε ( ∂ ∂ = − − ( ∂ ∂ ¸ ¸ ∂ ∂ ( = − − ( ∂ ∂ ¸ ¸ ∂ ( ∂ ∂ = − − ( ∂ ∂ ∂ ¸ ¸ (4.12) Dans ce cas, il n’est plus possible de séparer le système en deux sous-systèmes comme auparavant ; les cas TE z et TM z se trouvent donc mélangés et ne peuvent pas être traités séparément. IV.4 Réduction à une dimension Si nous considérons maintenant que ni le modèle ni les sources ne varient suivant les directions z et y, alors les dérivées spatiales des champs suivant ces deux directions sont nulles, et les deux systèmes se réduisent alors pour le mode transverse électrique TE zy 1 1 y z y y z z H E M t x H E J t x µ ε ∂ ∂ ( = − ( ∂ ∂ ¸ ¸ ∂ ( ∂ = − ( ∂ ∂ ¸ ¸ : (4.13) et pour le mode transverse magnétique TM zy 1 1 y z y y z z E H J t x E H M t x ε µ ∂ ∂ ( = − − ( ∂ ∂ ¸ ¸ ∂ ( ∂ = − − ( ∂ ∂ ¸ ¸ : (4.14) CHAPITRE IV : Méthode des différences finies dans le domaine temporel 79 IV.5 Algorithme de Yee L’idée principale de l’algorithme de la FDTD est de discrétiser les équations (4.4) et (4.5) dans leur forme différentielle et de les remplacer par un jeu d’équations aux différences finies. Il s’agit d’une méthode de numérisation qui permet de passer de l’expression analytique de l’équation à une approximation numérique. Cette méthode peut s’appliquer à toute dérivée partielle spatiale ou temporelle, du premier ou second ordre (développement en série de Taylor). IV.5.1 Principe des différences finies centrées D’un point de vue numérique, l’utilisation d’expressions programmables passe par la discrétisation des formulations analytiques. Les dérivées spatiales et temporelles des équations de Maxwell peuvent être traitées numériquement par la technique des différences finies. L’approximation des dérivées aux différents points de l’espace discret est réalisée par différenciation des valeurs des nœuds voisins ou point de dérivation. Figure 4.1 : Principe de calcul de la dérivée première de f(x) locale en x 0 . Soit f(x) une fonction continue représentant une composante du champ électrique ou magnétique et dérivable en point de l’espace comme le montre la figure 4.1. Les développements limités en série de Taylor à droite et à gauche de x 0 2 δ ± avec un décalage de s’écrivent : ( ) f x 0 ( / 2) f x δ + 0 ( ) f x 0 ( / 2) f x δ − 0 / 2 x δ − 0 / 2 x δ + 0 x x CHAPITRE IV : Méthode des différences finies dans le domaine temporel 80 2 3 0 0 0 0 0 1 1 ( ) ( ) ( ) ( ) ( ) . . . 2 2 2! 2 3 !2 f x f x f x f x f x δ δ δ δ | | | | ′ ′′ ′′′ + = + + + + | | \ . \ . (4.15) 2 3 0 0 0 0 0 1 1 ( ) ( ) ( ) ( ) ( ) . . . . 2 2 2! 2 3 !2 f x f x f x f x f x δ δ δ δ | | | | ′ ′′ ′′′ − = − + − + | | \ . \ . (4.16) En utilisant les équations (4.15) et (4.16), limitées à l’ordre 2, la dérivée première de f au point x 0 0 0 2 0 ( ) ( ) 2 2 ( ) ( ) f x f x f x O δ δ δ δ + − − ′ = + peut être évaluée de manière centrée à l’ordre 2 comme suit : (4.17) le terme 2 ( ) O δ représente l’erreur d’ordre 2 commise, et qui sera négligée par la suite. On appelle approximation centrée cette approximation de la dérivée. Les résultats qu’elle offre sont plus précis en comparaison avec ceux donnés par d’autres types d’approximations dites droits ou gauches, dont les formules (4.4) et (4.5) sont décrites respectivement ci-dessous : 0 0 0 ( ) ( ) 2 ( ) ( ) f x f x f x O δ δ δ + − ′ = + (4.18) 0 0 0 ( ) ( ) 2 ( ) ( ) f x f x f x O δ δ δ − − ′ = + (4.19) On remarque que le terme ( ) O δ qui est du premier ordre, est moins précis en comparaison avec celui du deuxième ordre de la dérivée centrée. Par conséquent, on utilisera l’approximation centrée dans notre étude pour discrétiser les dérivées partielles, spatiales et temporelles présentes dans les équations de Maxwell. En ajoutant les deux expressions (4.15) et (4.16) membre à membre, nous obtenons : 2 4 0 0 0 0 ( ) ( ) 2 ( ) ( ) ( ) 2 2 2 f x f x f x f x O δ δ δ δ | | ′′ + + − = + + | \ . (4.20) CHAPITRE IV : Méthode des différences finies dans le domaine temporel 81 où 4 ( ) O δ représente l’erreur de discrétisation, elle représente un terme négligeable devant les autres. Alors les opérateurs de dérivée première et deuxième peuvent être exprimés en différences finies selon : 0 0 0 ( ) ( ) 2 2 ( ) f x f x f x δ δ δ + − − ′ ≈ (4.21) 0 0 0 0 2 ( ) 2 ( ) ( ) 2 2 ( ) 2 f x f x f x f x δ δ δ + − + − ′′ ≈ | | | \ . (4.22) IV.5.2 Discrétisation des équations de Maxwell Une discrétisation spatiale et temporelle aux différences finies est effectuée pour la résolution des deux sous systèmes (4.4) et (4.5). La discrétisation des opérateurs de dérivation utilise un schéma centré des différences finies, avec une formulation dont l’erreur est du second ordre pour chaque pas de discrétisation (en d’autres termes : la diminution de moitié du pas de discrétisation réduit de 25% les erreurs d’évaluation des opérateurs de dérivation). a) La discrétisation spatiale : Pour montrer la manière dont on peut discrétiser les équations de Maxwell, reprenons l’une des six équations (4.4) et (4.5) avec M=0 : 1 y z x H E E t x z µ ∂ ∂ ∂ ( = − ( ∂ ∂ ∂ ¸ ¸ (4.23) Dans le calcul de H y , on fait intervenir la dérivée partielle de E x par rapport à z, donc d’après la définition de la dérivée centrée, le point où l’on calcule H y doit se trouver au milieu d’un segment parallèle à l’axe Oz, ayant comme extrémités deux points où E x est connu. De même, le calcul de H y fait intervenir la dérivée partielle de E z par rapport à x. Donc le point où l’on calcule H y se trouve également au milieu d’un segment parallèle à Ox ayant pour extrémités deux points où E z est connu. En conséquence, Hy représentée sur la figure 4.2 doit se trouver au milieu des deux points E x et des deux points E z . CHAPITRE IV : Méthode des différences finies dans le domaine temporel 82 Notons que les valeurs du champ électrique et du champ magnétique seront calculées en différents points du maillage. Elles seront respectivement appelées nœuds électriques et nœuds magnétiques. Figure 4.2 : Circulation du champ E autour de H Figure 4.3 : Disposition des nœuds électriques et des nœuds magnétiques dans le plan xOy On vient de voir que les nœuds H y doivent se trouver entre deux nœuds E x et entre deux nœuds E z ; de même les nœuds H x doivent se trouver entre deux nœuds E z et entre deux nœuds E y . Selon les équations (4.4) et (4.5), on constate également que le nœud E y , représenté sur la figure 4.4, doit se trouver entre deux nœuds H x et entre deux nœuds H z . Les nœuds E x doivent se trouver entre deux nœuds H y et entre deux nœuds H z . Les nœuds E z doivent se trouver entre deux nœuds H x et entre deux nœuds H y . CHAPITRE IV : Méthode des différences finies dans le domaine temporel 83 Figure 4.4 : Circulation du champ H autour de E L’arrangement des nœuds électriques et magnétiques doit donc respecter toutes ces conditions, et conduit au schéma de la maille de Yee représentée par la figure 4.5 suivante : Figure 4.5 : Cellule de YEE Les parallélépipèdes ou mailles élémentaires constituent le volume de calcul. Afin de le représenter selon le schéma décrit précédemment, on doit construire un maillage pour la structure étudiée. Précisons que dans le volume de calcul, sont toujours présents un nœud magnétique entre quatre nœuds électriques et un nœud électrique entre quatre nœuds magnétiques. Ainsi la dérivée centrée est utilisée pour toutes les dérivées spatiales présentes dans les équations de Maxwell. Pour représenter le volume de calcul, il est nécessaire de construire un maillage. CHAPITRE IV : Méthode des différences finies dans le domaine temporel 84 b) Construction du maillage Une discrétisation spatio-temporelle est nécessaire pour résoudre les équations de Maxwell (4.4) et (4.5). La discrétisation spatiale s’effectue dans un volume nécessairement fini. Dans le cas d’un maillage régulier, les dérivées spatiales sont évaluées dans les trois directions Ox, Oy, Oz avec des incréments constants : dx, dy, dz, appelés pas spatiaux. Ces derniers sont choisis par l’utilisateur et dépendent de la plus petite longueur d’onde présente dans la bande de fréquence d’analyse et de la géométrie de la structure à étudier. Le volume de calcul est donc un parallélépipède comme le montre la figure 4.6. Il est composé de (nx.ny.nz) cellules (ou mailles) élémentaires de taille dx , dy, dz. On va associer trois nœuds électriques et trois nœuds magnétiques pour chaque cellule élémentaire. Les valeurs du champ en ces nœuds seront notées Ex(i,j,k), Ey(i,j,k), Ez(i,j,k), Hx(i,j,k), Hy(i,j,k),Hz(i,j,k) et sont représentées ci-dessous dans la maille de Yee figure 4.6 où les entiers i, j, k représentent les indices de la cellule dans le maillage et varient respectivement de 1 à nx, 1 à ny, et 1 à nz. Figure 4.6 : Extraction d’une cellule élémentaire c) Discrétisation temporelle : Considérons de nouveau la même équation (4.23). Elle fait apparaître dans le membre de gauche la dérivée du champ magnétique par rapport au temps, tandis que le membre de droite est considéré à un instant t. Si on prend en considération le principe de la dérivée centrée, on en déduit que le membre de droite (le champ électrique) doit être calculé entre deux instants successifs où on calcule le membre de gauche (le champ magnétique). En tenant compte des six équations de Maxwell, on arrive à la conclusion que le champ électrique et le champ magnétique ne doivent pas être calculés aux mêmes instants, mais à des instants décalés. CHAPITRE IV : Méthode des différences finies dans le domaine temporel 85 Pour le cas d’une discrétisation temporelle uniforme, avec un pas d’échantillonnage dt, le champ électrique sera calculé pour des multiples impairs de dt/2, et le champ magnétique pour les multiples pairs de dt/2 comme le montre la figure 4.7 : Figure 4.7 : Calcul de H à l’instant ndt et calcul de E à l’instant (n+0.5)dt IV.5.3 Equations de Maxwell aux différences centrées En utilisant le principe des différences finies centrées et les notations de Kane Yee, pour une fonction , , ( , , , ) n i j k u i x j y k z n t u ∆ ∆ ∆ ∆ = où i, j, k, n sont des entiers, la dérivée temporelle de u au point (i, j, k) s’exprime alors simplement : 1/ 2 1/ 2 , , , , 2 ( , , , ) ( ) n n i j k i j k u u u i x j y k z n t O t t t + − − ∂ ( ∆ ∆ ∆ ∆ = + ∆ ¸ ¸ ∂ ∆ (4.24) de même que ses dérivées spatiales à l’instant nΔt : 1/ 2, , 1/ 2, , 2 , 1/ 2, , 1/ 2, 2 , , 1/ 2 , , 1/ 2 2 ( , , , ) ( ) ( , , , ) ( ) ( , , , ) ( ) n n i j k i j k n n i j k i j k n n i j k i j k u u u i x j y k z n t O x x x u u u i x j y k z n t O y y y u u u i x j y k z n t O z z z + − + − + − − ∂ ( ∆ ∆ ∆ ∆ = + ∆ ¸ ¸ ∂ ∆ − ∂ ( ∆ ∆ ∆ ∆ = + ∆ ¸ ¸ ∂ ∆ − ∂ ( ∆ ∆ ∆ ∆ = + ∆ ¸ ¸ ∂ ∆ (2.25) En appliquant le principe des différences finies centrées avec ces notations, une approximation numérique du système des équations différentielles couplées (4,4) et (4,5) peut maintenant être obtenue: , 1, 1/ 2 , , 1/ 2 1/ 2 1/ 2 , 1/ 2, 1/ 2 , 1/ 2, 1/ 2 , 1/ 2, 1 , 1/ 2, , 1/ 2, 1/ 2 , 1/ 2, 1/ 2 , 1/ 2, 1/ 2 1 n n z z i j k i j k n n n n y y x x i j k i j k i j k i j k i j k n i j k x i j k H H y H H E E t z E ε σ + + + + − + + + + + + + + + + + + + ( − ( ∆ ( ( − ( − = − ( ∆ ∆ ( ( − ( ( ( ¸ ¸ (4.26) CHAPITRE IV : Méthode des différences finies dans le domaine temporel 86 1/ 2, 1, 1 1/ 2, 1, 1/ 2 1/ 2 1/ 2, 1, 1/ 2 1/ 2, 1, 1/ 2 , 1, 1/ 2 1, 1, 1/ 2 1/ 2, 1, 1/ 2 1/ 2, 1, 1/ 2 1/ 2, 1, 1/ 2 1 n n x x i j k i j k n n n n y y z z i j k i j k i j k i j k i j k n i j k y i j k H H z E E H H t x E ε σ − + + − + + − − + + − + + + + − + + − + + − + + − + + ( − ( ∆ ( ( − − ( = − ( ∆ ∆ ( ( − ( ( ¸ ¸ (4.27) , 1/ 2, 1 1, 1/ 2, 1 1/ 2 1/ 2 1/ 2, 1/ 2, 1 1/ 2, 1/ 2, 1 1/ 2, 1, 1 1/ 2, , 1 1/ 2, 1/ 2, 1 1/ 2, 1/ 2, 1 1/ 2, 1/ 2, 1 1 n n y y i j k i j k n n n n z z x x i j k i j k i j k i j k i j k n i j k z i j k H H x E E H H t y E ε σ + + − + + + − − + + − + + − + + − + − + + − + + − + + ( − ( ( ∆ ( − − ( = − ( ∆ ∆ ( ( − ( ( ¸ ¸ (4.28) 1/ 2 1/ 2 1/ 2, 1, 3/ 2 1/ 2, 1, 1/ 2 1 1/ 2 1/ 2 1/ 2, 1, 1 1/ 2, 1, 1 1/ 2, 1, 1 1/ 2, 1/ 2, 1 1/ 2, 1, 1 1/ 2 * 1/ 2, 1, 1 1/ 2, 1, 1 1 n n y y i j k i j k n n n n x x z z i j k i j k i j k i j k i j k n i j k x i j k E E z H H E E t y H ε σ + + − + + − + + + + + − + + − + + − + + − + + − + + + − + + − + + ( − ( ( ∆ ( − − ( = − ( ∆ ∆ − ¸ ¸ ( ( ( ( (4.29) 1/ 2 1/ 2 1/ 2, 1/ 2, 1 1/ 2, 1/ 2, 1 1 1/ 2 1/ 2 , 1/ 2, 1 , 1/ 2, 1 , 1/ 2, 3/ 2 , 1/ 2, 1/ 2 , 1/ 2, 1 1/ 2 * , 1/ 2, 1 , 1/ 2, 1 1 n n z z i j k i j k n n n n y y x x i j k i j k i j k i j k i j k n i j k y i j k E E x H H E E t z H ε σ + + + + + − + + + + + + + + + + + + + + + + + + + + ( − ( ∆ ( ( − − ( = − ( ∆ ∆ ( ( − ( ( ¸ ¸ (4.30) CHAPITRE IV : Méthode des différences finies dans le domaine temporel 87 1/ 2 1/ 2 , 3/ 2, 1/ 2 , 1/ 2, 1/ 2 1/ 2 1/ 2 1 , 1, 1/ 2 , 1, 1/ 2 1/ 2, 1, 1/ 2 1/ 2, 1, 1/ 2 , 1, 1/ 2 1/ 2 * , 1, 1/ 2 , 1, 1/ 2 1 n n x x i j k i j k n n n n y y z z i j k i j k i j k i j k i j k n i j k z i j k E E y E E H H t x H ε σ + + + + + + + + + + + + + + + + − + + + + + + + + + ( − ( ∆ ( ( − ( − = − ( ∆ ∆ ( ( − ( ( ( ¸ ¸ (4.31) IV.5.4 Dispersion numérique Le choix fixé des paramètres de discrétisation en différences finies induit des erreurs lors de la propagation (dispersion non physique des signaux qui se propagent sur la grille de calcul). En effet, le spectre de la source dont l’amplitude est notable sur un intervalle borné, n’est qu’en partie correctement discrétisé : l’idée étant de définir les pas spatiaux et temporels pour la fréquence la plus énergétique dite fréquence de pic de l’impulsion. Au delà d’une certaine fréquence, chaque maille agit par conséquent comme un réflecteur. Ce phénomène apparaît en général lorsque la longueur d’onde associée à la fréquence correspond à moins de 3 mailles, et la grille de Yee se comporte par conséquent comme un passe bas. Au bout d’un certain temps, on observe ainsi un train d’onde haute fréquence en retard par rapport à l’ensemble (Fig. 4.8). Figure 4.8 : Dispersion numérique pour une impulsion de type gaussienne; on remarque la distorsion en amont de l'impulsion gaussienne La mise en évidence théorique de la dispersion numérique se fait en comparant l’expression discrétisée et l’expression analytique de l’expression entre le nombre d’onde et la pulsation, dite relation de dispersion. Pour une grille 3D à maillage orthogonal, et dans le cas d’une onde plane incidente avec un angle α, on a la relation de dispersion suivante : CHAPITRE IV : Méthode des différences finies dans le domaine temporel 88 2 2 2 2 1 1 1 1 sin( ) sin( ) sin( ) sin( ) 2 2 2 2 y x z k y k x k z t t x y z ω ∆ ( ∆ ∆ ∆ ( ( ( = + + ( ( ( ( ∆ ∆ ∆ ∆ ¸ ¸ ¸ ¸ ¸ ¸ ¸ ¸ (4.32) L’écriture en coordonnées sphériques de cette équation et pour un pas de discrétisation spatial isotrope Δ= Δx = Δy = Δz, donne : 2 2 2 2 2 cos( ) sin( ) sin( ) sin sin sin 2 2 2 sin( ) cos( ) sin 2 t k k t k ω ϕ θ ϕ ϕ θ ∆ ∆ ∆ ∆ | | | | | | = + | | | ∆ \ . \ . \ . ∆ | | + | \ . (4.33) La résolution de l’équation (4.33) dans laquelle k est l’inconnue s’effectuée de manière itérative par la méthode de Newton : 2 2 2 1 sin ( ) sin ( ) sin ( ) sin(2 ) in(2 ) sin(2 ) i i i i i i i i Ak Bk Ck D k k A Ak B Bk C Ck + + + − = − + + (4.34) avec 2 cos( ) sin( ) sin( ) sin( ) cos( ) , B= , C= , D= sin 2 2 2 2 t A t ϕ θ ϕ ϕ θ ω ∆ ∆ ∆ ∆ ∆ | | = | ∆ \ . (4.35) La longueur d’onde dans le milieu discrétisé est 2 / g k λ π = . En choisissant / 2 t c ∆ = ∆ , on obtient une mesure de la dispersion numérique dans le rapport 0 p g v c λ λ = où 0 λ est la longueur d’onde dans le vide, v p est la vitesse de phase dans le milieu discrétisé et c la vitesse dans le vide. Pour un pas de discrétisation spatiale égal à un dixième de la longueur d’onde, on observe une dispersion numérique inférieure à 2%. Cependant la dispersion numérique augmente linéairement en fonction de la distance parcourue par l’onde. Par exemple, une onde plane se propageant sur une distance approchant 10 fois sa longueur d’onde, cumule environ 10% d’erreur sur la vitesse de phase [5]. Afin de limiter la dispersion numérique, l’augmentation de l’ordre du schéma de discrétisation FDTD pourrait être une solution efficace, mais aurait un coût de calcul non négligeable. CHAPITRE IV : Méthode des différences finies dans le domaine temporel 89 IV.5.5 Critères de convergence et de stabilité de l’algorithme La convergence du schéma numérique est assurée si la vitesse de propagation numérique de l’onde dans la grille est finie et supérieure ou égale à la vitesse de phase de l’onde réelle. L’application de cette contrainte implique une relation entre les pas de discrétisation temporelle et spatiale donnés par la condition de Courant-Friedrichs-Lewy [6] qui s’écrit (à une dimension) : max Num P h V V t = ≥ ∆ Soit max 1 P V t h ∆ ≤ (4.36) où h représente l’incrément spatial et Δ t l’incrément temporel, V Pmax ( ) ( ) ( ) max 2 2 2 1 1 1 1 P t V x y z ∆ ≤ + + ∆ ∆ ∆ représente la vitesse de phase maximale de l’onde dans le volume de calcul. Dans le cas général, pour un maillage cartésien en 3D, cette condition devient : (4.37) et pour un maillage cubique où δ = Δx = Δy = Δz, l’équation (4.37) se réduit à : max 3 P t V δ ∆ ≤ (4.38) Dans ce cas, la condition sur le nombre d’étapes nécessaires pour une période d’observation T, telle que : T n t = ∆ , s’écrit : max 3 P V T n δ ≥ (4.39) Cette équation montre que pour les applications à basse fréquence, l’importance de la période d’observation implique un nombre plus important d’itérations. La précision des résultats est principalement liée à la finesse des pas spatio-temporels. Par ailleurs, pour obtenir un régime établi, le calcul doit se faire sur un minimum de cinq périodes du signal source et les données correspondantes doivent être enregistrées. Ces constatations impliquent une capacité mémoire importante et un temps de calcul élevé. CHAPITRE IV : Méthode des différences finies dans le domaine temporel 90 IV.6 Sources et signaux d'excitation Le choix de l’excitation dépend de plusieurs facteurs parmi lesquels, le type de la structure à étudier et la bande de fréquences ciblée. L’excitation se traduira par un signal numérique qui va se propager dans la structure. Cette variation imposée à un endroit approprié du maillage, a une forme, une durée et un emplacement particulier. Numériquement, on peut choisir une forme arbitraire pour l’excitation. On a toutefois intérêt pour des problèmes de convergence des résultats, à choisir une excitation proche de la forme du champ réel dans la structure. IV.6.1 Impulsion Gaussienne Il est souhaitable d’utiliser une excitation capable de remplir certaines conditions, comme une étude sur une large bande spectrale, d’une durée temporelle raisonnable, continue, et facilement interprétable. A partir de toutes ces données, on peut dire que l’excitation la mieux adaptée est une Gaussienne dont l’équivalent fréquentiel est une « demi-gaussienne ». En effet, son expression analytique est simple; le spectre en fréquence est facilement contrôlable. Le signal est borné dans le temps, son évolution est lisse et ne présente pas de variations trop rapides qui pourraient générer des erreurs de calcul. Le fait d’utiliser la Gaussienne permet en une simulation de connaître, en faisant une transformée de Fourier, la réponse sur une large bande de fréquence. Il est difficile, dans la réalité, de reproduire une telle excitation de façon expérimentale. Mais la forme choisie n’intervient que pour la simulation, comme un intermédiaire de calcul qui permet de connaître la réponse du système sur une large bande de fréquence. Une source gaussienne sera définie de la façon suivante : 2 0 0 2 ( ) ( ) exp t t e t A T | | − = − | \ . (4.40) où A 0 t n t = ∆ représente l’amplitude de la gaussienne, (n est le nombre d’itérations et t ∆ est le pas temporel), T 0 désigne le retard par rapport à l’instant t = 0, T est proportionnelle à la largeur à mi-hauteur de la gaussienne. T et la fréquence maximale de la bande étudiée f max max 1 2 T f ≈ sont reliées par l’expression: (4.41) CHAPITRE IV : Méthode des différences finies dans le domaine temporel 91 Figure 4.9 : Allure temporelle de la gaussienne La dérivée gaussienne (Fig. 4.9) est aussi une forme d’excitation connue, son expression temporelle e(t) se définie par : 2 0 0 0 2 ( ) ( ) 2 ( ) exp t t t t e t A T T | | − − = − − | \ . (4.42) Figure 4.10 : Allure temporelle d’une dérivée gaussienne IV.6.2 Excitation sino-gaussienne L’excitation sino-gaussienne (gaussienne modulée par un sinus) permet une modélisation du continu jusqu'à une fréquence maximale. Il peut s'avérer nécessaire de modéliser une bande de fréquences n'incluant pas le continu (Guides d’ondes en bande X…). Pour ce faire, il suffit de multiplier la gaussienne par une sinusoïde dont la fréquence va correspondre à la fréquence centrale de la bande spectrale à étudier. 2 0 0 2 ( ) ( ) sin(2 ) exp n t T S n f n t T π | | ∆ − = ∆ − | \ . (4.43) CHAPITRE IV : Méthode des différences finies dans le domaine temporel 92 f 0 Figure 4.11 : Forme temporelle d'une source sinusoïdale modulée par une gaussienne IV.6.3 Excitation par une onde plane représente la fréquence centrale de la bande étudiée. La largeur de la bande de fréquence étudiée est environ égale à 1/T. L’utilisation d’une onde plane est très courante dans les modélisations électromagnétiques et son implémentation dans les logiciels est indispensable au traitement des problèmes liés à l’illumination de structures. La méthode FDTD utilise une formulation champ total/champ diffracté [7] qui se base sur l’utilisation du concept des surfaces de d’Huygens [8]. Le domaine de calcul est divisé en deux parties (Figure 4.12). Les champs incidents de l’onde plane sont introduits sur une surface virtuelle séparant les deux zones (TF/SF), de manière à ce qu’ils soient confinés dans la zone de champ total. Dans ce dernier, la FDTD prend en compte la somme du champ incident de l’onde plane et du champ diffracté par les objets, alors que seul ce dernier est propagé dans la zone de champ diffracté. Le champ source tient compte de la lumière diffractée dans le reste de la zone étudiée. Cela évite l’apparition de réflexions non-physiques au niveau de la source. La zone associée au champ diffracté seule n’a pas de réalité physique. Pour rendre opérationnelle ce type de source, des corrections du champ calculé par l’algorithme de Yee sont nécessaires au niveau des frontières entre le champ total et le champ diffracté. On appelle « champ total » la somme du champ diffracté par la structure et du champ incident dû à la source. Dans ce cas, la source est introduite à la transition entre zone de CHAPITRE IV : Méthode des différences finies dans le domaine temporel 93 champ total et de champ diffracté. La séparation du volume de calcul en deux zones distinctes utilise la propriété de linéarité des équations de Maxwell. t s d E E E = + (4.44) t s d H H H = + (4.45) Figure 4.12 : Introduction de la source entre les deux zones du TF/SF (cas 2D). avec s E et s H représentant les valeurs des champs incident supposés connus, d E et d H les valeurs des champs diffusés initialement inconnus. IV.7 Conditions d’absorption aux limites Les équations de Maxwell sont résolues dans un domaine de calcul dont les dimensions sont nécessairement finies. Néanmoins certaines simulations numériques demandent des conditions d’espace libre. Il faut donc soit agrandir le domaine de calcul de telle sorte que les ondes réfléchissantes ne perturbent pas les résultats, soit appliquer des conditions particulières sur les frontières afin d’obtenir un domaine non borné. La première solution est restrictive, elle va demander beaucoup trop de place mémoire défavorisant ainsi le temps de calcul, par contre la deuxième est la plus avantageuse pour les simulations numériques. CHAPITRE IV : Méthode des différences finies dans le domaine temporel 94 IV.7.1 Bref état de l’art Pour absorber les ondes sortantes, diverses méthodes ont été utilisées dans les codes numériques dont la première qui fut la méthode RBC (Radiation Boundary Conditions). C’est une famille d’opérateurs aux dérivés partielles en temps et espace basée sur l’équation d’Helmoltz. B. Engquist et A. Majda [9] ont proposé une solution qui reprend l’opérateur de A. Bayliss et E. Turkel [10] en l’adaptant aux coordonnées cartésiennes. L’équation d’Helmoltz est ainsi divisée suivant chaque axe en deux opérateurs correspondant aux deux sens de propagation. Puis G. Mur [11] et a étendu l’opérateur pour qu’il devienne un opérateur de second ordre. La couche adaptée est une méthode qui consiste à entourer le domaine de calcul par un milieu absorbant dont l’impédance est égale à celle du vide. Les travaux de R. Holland [12] présentent cette approche et cette nouvelle famille appelée ABC (Absorbing Boundary Conditions). Néanmoins, aucune de ces techniques d’absorption d’ondes réfléchissantes n’est parfaite. Elles ne sont parfaitement efficaces que pour des cas particuliers. Seules les ondes planes se propageant en incidence normale seront absorbées. De plus, ces techniques nous imposent une contrainte géométrique bien connue comme la nécessité de placer les frontières à une certaine distance de la structure à étudier (de l’ordre de min 3λ ). Mais une nouvelle technique de simulation d’espace libre va considérablement accélérer l’efficacité des ABCs. J-P Bérenger va, en effet, créer en 1994 [13] les PMLs (Perfectly Matched Layers). Cette technique est basée sur l’utilisation de couches absorbantes de R. Holland dont le milieu est remplacé par un nouveau milieu spécialement élaboré pour absorber sans réflexion les ondes électromagnétiques quelque soit l’angle d’incidence. Ce nouveau milieu est composé de couches parfaitement adaptées quelque soit l’angle d’incidence et quelque soit la fréquence. De plus, la structure à étudier peut très bien être placée proche des PMLs sans qu’il y ait de perturbations trop fortes. Les résultats d’absorptions obtenus par les codes numériques utilisant les PMLs sont bien meilleurs que ceux obtenus par les autres méthodes. Néanmoins, J-P. Bérenger présente en 1999 la théorie des ondes évanescentes dans les milieux homogènes et discrétisés. Il montre qu’il existe de fortes réflexions numériques dues CHAPITRE IV : Méthode des différences finies dans le domaine temporel 95 aux ondes évanescentes lorsqu’une structure (de type plaque) se trouve très proche des PMLs. Pour pallier à ce problème, les PMLs sont réécrites en se basant sur les « Stretched Coordinate Space » [14]; c’est une nouvelle formulation des équations de Maxwell où le repère cartésien est étendu en coordonnée complexe. Puis M. Kuzuoglu et R. Mittra [15] proposent une forme de PML strictement causal, continu en décomposant l’opérateur en une partie réelle et une partie imaginaire. Nous faisons référence ici au CFS-PML (Complex Frequency Shifted-PML). C’est une nouvelle forme de PML qui sera construite permettant d’optimiser l’efficacité de l’absorption des PMLs classiques. Les travaux montrent qu’elles sont très performantes pour les ondes évanescentes. Une technique d’implantation dans un code FDTD est introduite par J. Alan Roden et Stephen D. Gedney [16] basée sur la formulation des « Stretched Coordinate Space » et la convolution récursive. Cette nouvelle adaptation de PML est appliquée dans la méthode des Différences Finies et elle est référencée par CPML (Convolution Perfectly Matched Layers). Il est montré que cette technique est robuste et fiable. Les CPMLs absorbent toutes les ondes électromagnétiques aussi bien dans les milieux stratifiés, à pertes ou dispersifs. IV.7.2 Conditions périodiques aux limites Les conditions périodiques aux limites (Periodic Boundary Conditions, PBCs) sont un artifice de simulation qui permet de s’affranchir des effets de bords liés à la troncature de l’espace de calcul. Elles sont principalement utilisées dans le domaine de la chimie, mais on peut les employer pour les problèmes de propagation d’ondes électromagnétiques. Le principe consiste à considérer que le domaine de calcul est répété une infinité de fois dans l’espace, apposé à lui-même, et donc que des relations particulières existent entre les valeurs traitées sur les différentes faces du domaine de base. Afin de réussir un pavage parfait de l’espace, il importe que le domaine élémentaire de calcul soit de forme adéquate. En général, on s’oriente vers le parallélépipède ou le prisme rectangulaire. Ce sont les choix les plus intuitifs, mais ils sont toutefois plus coûteux en termes de maillage spatial que l’utilisation de l’octaèdre tronqué comme forme de base. Dans le cas d’un espace parallélépipédique, les faces du pavé seront liées entre elles deux à deux, car chacune sera, dans la répétition du motif, adjointe à la face opposée ; il s’agira au final d’une seule et même interface. Si l’on considère l’évolution d’un champ électromagnétique, les composantes du champ aux frontières subiront le même type de CHAPITRE IV : Méthode des différences finies dans le domaine temporel 96 contraintes. Il n’y aura pas de réflexion parasite à l’interface, car toute onde incidente sortant du domaine de calcul se verra confrontée, par périodicité spatiale, au même type de milieu que celui d’où elle provient – ceci supposant qu’une certaine adéquation de caractéristiques électromagnétiques existe entre les milieux bordant les faces opposées ; l’intensité de la réflexion parasite à l’interface sera donc fonction de cette conformité. Si l’adaptation d’impédance peut être parfaite avec la méthode de la CPL, il n’en subsiste pas moins quelques problèmes, surtout pour de petits espaces de calcul. En effet, la répétition du motif est seulement virtuelle, et finalement on ne travaille que sur une seule entité spatiale. En d’autres termes, si une onde incidente quitte le domaine de calcul par la frontière gauche, elle réapparaîtra à la frontière droite. On voit immédiatement la limite du modèle : il peut se révéler judicieux lorsque l’on dispose, par exemple, de flux continus de particules, tels que rencontrés en chimie. En revanche, pour la simulation d’une propagation d’onde électromagnétique, il est non physique de modéliser la transposition d’un phénomène d’un point à l’autre de l’espace instantanément. Le seul contexte dans lequel cette approximation se montrerait valable, serait l’utilisation d’un espace de calcul très grand devant la zone étudiée, ce qui permettrait d’occulter le retour d’ondes transposées. Mais dans de telles conditions, on peut également s’interroger sur l’opportunité de définir des conditions aux limites. IV.7.3 Conditions d’Engquist-Majda-Mur Les conditions aux limites développées par Mur [11] se fondent sur des résultats et équations énoncés par Engquist et Majda en 1977 [9]. La condition n’est opérationnelle que dans le cas d’un maillage cartésien et le principe repose sur la factorisation des opérateurs de dérivées partielles dans l’équation d’onde. Ces conditions aux limites permettent d’obtenir une réflexion aux frontières de l’ordre de 1/100 (rapport des amplitudes). L’équation d’onde s’écrit pour une onde électromagnétique sous la forme : 2 2 2 2 2 2 2 2 2 1 0 u u u u x y z c t ∂ ∂ ∂ ∂ + + − = ∂ ∂ ∂ ∂ (4.46) où u est une composante quelconque parmi les six composantes du champ électromagnétique. L’équation peut se récrire en utilisant un opérateur L appliqué à la composante u, avec : 0 Lu = et 2 2 2 2 2 2 2 2 2 1 u u u u L x y z c t ∂ ∂ ∂ ∂ = + + − ∂ ∂ ∂ ∂ (4.47) CHAPITRE IV : Méthode des différences finies dans le domaine temporel 97 On peut écrire l’opérateur L sous la forme du produit de deux opérateurs x L + et x L − , où : 2 2 1 1 1 1 y z x c t c t L x c t ∂ ∂ ∂ ± ∂ ∂ ∂ ∂ ∂ | | | | ∂ ∂ = ± − − | | ∂ ∂ \ . \ . (4.48) Cette factorisation induit la considération d’une propagation le long de l’axe des x (d’où l’indice x). Engquist et Majda ont montré en que l’égalité fonctionnelle 0 L u ± = déterminait les conditions sur les composantes tangentielles du champ pour éliminer l’onde réfléchie en bordure du domaine de calcul. Considérons une onde plane se propageant dans le sens des x croissants et incidente sur la paroi x = d. Autour de l’incidence normale, on peut considérer que les dérivées suivant y et z sont petites par rapport à la dérivée temporelle, et donc en notant : 1 1 , S y z y z c t c t S ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ = = (4.49) on peut écrire : 2 2 2 2 1 1 1 1 2 2 y z y z S S S S − − ≈ − − (4.50) L’opérateur L s’écrit alors : 2 2 1 1 1 1 2 2 x y z L S S x c t + ∂ ∂ | | = ± − − | ∂ ∂ \ . (4.51) et la condition 0 L u ± = : 2 2 2 2 2 2 2 1 2 2 x u u c u c u L u x t c t y z + ∂ ∂ ∂ ∂ = + − − ∂ ∂ ∂ ∂ ∂ (4.52) Cette équation exprime les conditions de Mur au second ordre. À l’examen de cette équation, on peut formuler deux remarques : • La solution n’est strictement valable que pour une incidence normale. Des réflexions parasites apparaîtront lorsqu’on s’en écartera. • Les conditions sont les mêmes pour les champs électrique et magnétique. CHAPITRE IV : Méthode des différences finies dans le domaine temporel 98 IV.7.4 Couches parfaitement adaptées « PML » Ces conditions aux limites sont certainement les conditions d’absorption les plus performantes aujourd’hui. Elles permettent de descendre à des réflexions en amplitude de l’ordre de 10 -5 σ , sur une très large gamme d’incidences et de fréquences. La technique PML (Perfectly Matched Layers) repose sur le principe d’adaptation d’impédance à l’interface entre deux milieux de même indice mais dont l’un est absorbant avec des conductivités électrique et magnétiques * σ non nulles (Fig. 4.13). Cette condition d’adaptation s’écrit : * 0 σ σ ε µ = (4.53) Figure 4.13 : Principe d’adaptation d’impédance L’onde arrivant du milieu incident n’est pas réfléchie vers celui-ci et se trouve atténuée dans le milieu absorbant. Mais dans ce cas, l’adaptation d’impédance n’est possible qu’à incidence normale, des réflexions parasites à l’interface apparaissent dans le cas où l’onde arrive à incidence oblique. Pour y remédier, Berenger a proposé un milieu absorbant artificiellement biaxe. L’absorption est non nulle suivant la normale à l’interface entre les deux milieux et elle est nulle suivant l’axe parallèle à l’interface (Figure 4.14). Dans le milieu PML, l’onde plane incidente est décomposée fictivement en deux ondes : 1) Une onde se propageant sous incidence normale et vérifie l’équation 4.53. elle est atténuée et absorbée par le milieu PML et ne subit qu’une très faible réflexion vers le milieu incident. CHAPITRE IV : Méthode des différences finies dans le domaine temporel 99 2) Une deuxième onde à incidence rasante qui ne subit aucune absorption dans le milieu PML. Cette onde, se propageant parallèlement à l’interface entre les deux milieux ne subit aucune réflexion et voit un milieu identique à celui de la fenêtre principale. Figure 4.14 : Principe de fonctionnement d’une PML La forte absorption de l’onde dans le milieu PML, peut engendrer des réflexions parasites vers la structure étudiée. Ces réflexions, qui sont d’origine purement numérique, proviennent de la discontinuité induite par la discrétisation spatiale. En effet, la technique de discrétisation aux différences centrées est inadaptée dans le cas où les champs subissent de fortes variations. Une solution à ce problème est d’imposer une augmentation progressive de façon polynomiale de l’absorption σ dans la couche PML. Elle est donnée comme suit : max m pml x e σ σ | | = | \ . (4.54) avec max σ est la conductivité maximale, pml x représente la profondeur dans la région PML mesurée à partir de l’interface, e désigne l’épaisseur de la couche PML et m dénote le degré de la loi polynomiale qui est généralement égal 2. L’absorption part de zéro à l’interface PML-domaine de calcul et augmente pour atteindre sa valeur maximale au bord extérieur de la PML. La condition de mur électrique qui consiste à forcer à zéro les composantes tangentielles du champ électrique est appliquée sur le bord extérieur de la couche PML. Dans le milieu PML, les composantes électriques et magnétiques sont dédoublées. Au total, on obtient 12 composantes électromagnétiques qui sont décrites par [5, 17, 18] : CHAPITRE IV : Méthode des différences finies dans le domaine temporel 100 ( ) y xy zx zy E H H t y ε σ ∂ ∂ | | + = + | ∂ ∂ \ . (4.55a) ( ) z xz yx yz E H H t z ε σ ∂ ∂ | | + = − + | ∂ ∂ \ . (4.55b) ( ) z yz xy xz E H H t z ε σ ∂ ∂ | | + = + | ∂ ∂ \ . (4.55c) ( ) x yx zx zy E H H t x ε σ ∂ ∂ | | + = − + | ∂ ∂ \ . (4.55d) ( ) x zx yx yz E H H t x ε σ ∂ ∂ | | + = + | ∂ ∂ \ . (4.55e) ( ) y zy xy xz E H H t z ε σ ∂ ∂ | | + = − + | ∂ ∂ \ . (4.55f) * ( ) y xy zx zy H E E t y µ σ ∂ ∂ | | + = − + | ∂ ∂ \ . (4.56a) * ( ) z xz yx yz H E E t z µ σ ∂ ∂ | | + = + | ∂ ∂ \ . (4.56b) * ( ) z yz xy xz H E E t z µ σ ∂ ∂ | | + = − + | ∂ ∂ \ . (4.56c) * ( ) x yx zx zy H E E t x µ σ ∂ ∂ | | + = + | ∂ ∂ \ . (4.56d) * ( ) x zx yx yz H E E t x µ σ ∂ ∂ | | + = − + | ∂ ∂ \ . (4.56e) * ( ) y zy xy xz H E E t y µ σ ∂ ∂ | | + = + | ∂ ∂ \ . (4.56f) Pour résoudre les équations précédentes dans le milieu PML, il faut les discrétiser aux différences centrées dans l’espace et dans le temps. CHAPITRE IV : Méthode des différences finies dans le domaine temporel 101 Figure 4.15 : Domaine de calcul FDTD entouré de PML Notons aussi que les PML de type Bérenger ont deux limitations importantes : d’une part, elles n’absorbent pas les ondes évanescentes, et d’autre part, elles ne sont pas adaptées à la simulation de milieux dispersifs. D’autres modèles de conditions aux limites comme celles de type UPML (Uniaxial Perfectly Matched Layer) ou CPML (Convolution Perfectly Matched Layer) ayant une interprétation physique ont été développées comme ceux qui utilisent des matériaux anisotropes et qui permettent d’absorber les ondes évanescentes et qui simulent des milieux dispersifs. IV.8 Implémentation des milieux dispersifs L’implémentation des milieux dispersifs, décrits soit par le modèle de Debye, soit par le modèle de Lorentz, peut se faire essentiellement de deux manières distinctes : par convolution récursive (Recursive Convolution method) ou par Transformée de Fourier (Auxiliary Differential Equation method) [5]. IV.8.1 Méthode RC Luebbers et al [19] intègrent un modèle de dispersion linéaire par convolution récursive RC (Recursive Convolution) dans la formulation par différences finies. Soit la relation dans le domaine spectral entre le champ électrique E , le flux correspondant D et la permittivité diélectrique ε : CHAPITRE IV : Méthode des différences finies dans le domaine temporel 102 ( ) ( ) ( ) D E ω ε ω ω = (4.57) En décomposant la permittivité en un terme statique et en un terme de polarisation on obtient : | | 0 0 ( ) ( ) ( ) ( ) ( ) ( , ) ( ) e D r r E r P r r r E r ε ε ε ε χ ω ∞ ∞ = + = + (4.58) avec ε ∞ la permittivité haute fréquence. Dans le cadre d’une susceptibilité électrique χ e 0 0 ( , ) ( ) ( , ) ( , ) * ( , ) e D r t r E r t r t E r t ε ε ε χ ∞ = + caractérisant une dispersion linéaire (linéarité concernant la relation entre la polarisation P et E), la transformation de Fourier inverse donne : (4.59) avec * l’opérateur de convolution temporelle. En insérant (4.59) dans (4.2), on obtient : 0 0 ( , ) ( , ) ( , ) * ( , ) ( ) ( , ) e r t E r t H r t E r t r E r t t t χ ε ε ε σ ∞ ∂ ∂ ∇× = + + ∂ ∂ (4.60) Pour un modèle de type Debye, l’expression de la susceptibilité est de la forme suivante : 0 ( , ) 1 s e r j t ε ε χ ω ω ∞ − = + (4.61) où s ε est la permittivité relative basse fréquence, ε ∞ 0 t est la permittivité haute fréquence et est un temps de relaxation caractéristique du milieu. Soit, ( , ) ( ) t e r t Ge U t γ χ − = (4.62) U(t) est la fonction unité, 0 ( ) / s G t ε ε ∞ = − et 0 1/ t γ = Ce qui donne d’après (4.60) : 0 0 ( , ) ( , ) ( ) * ( , ) ( ) ( , ) n E r t H r t r E r t r E r t t ε ε ε ψ σ ∞ ∂ ∇× = + + ∂ (4.63) CHAPITRE IV : Méthode des différences finies dans le domaine temporel 103 où la fonction ψ à l’instant n est définie telle que : 1 0 0 ( ) exp( ( ) ) ( ) n m n m r G n m t E r ψ ε γ γ − = = − − ∆ ∑ (4.64) On en déduit la relation de récurrence obtenue grâce aux propriétés de la fonction exponentielle : | | 1 ( ) ( ) exp( ) 1 exp( ) ( ) n n n r r t G t E r ψ ψ γ γ + = − ∆ + − − ∆ (4.65) Cette relation permet d’éviter le stockage de toutes les valeurs de champ pendant les itérations, puisque la relation de convolution est évaluée dans son intégralité à chaque pas de temps. IV.8.2 Méthode ADE La méthode ADE (Auxiliary Differential Equation) [5, 20] est basée sur les propriétés d’inversion de la transformée de Fourier entre les domaines spectral et temporel. Ces mêmes propriétés permettent d’établir les équivalences suivantes : 1 1 2 n n E E E j E t t ω + − ∂ − ⇔ = ∂ ∆ (4.66) ( ) 2 1 1 2 2 2 2 n n n E E E E E t t ω + − ∂ − + − ⇔ = ∂ ∆ (4.67) Considérons à titre d’illustration le modèle de dispersion sur la permittivité de type Debye : 0 0 1 s j t ε ε ε ε ε ω ∞ ∞ | | − = + | + \ . (4.68) En replaçant (4.68) dans (4.57), un peu d’algèbre suivi d’une formulation aux différences finies permet d’obtenir pour la composante électrique suivant x : 1 1 1 2 3 n n n n x x x x E a D a D a E + + = + + (4.69) avec 0 1 0 2 2 s t t a t t ε ε ∞ ∆ + = + ∆ , 0 2 0 2 2 s t t a t t ε ε ∞ ∆ − = + ∆ et 3 0 2 2 s s t t a t t ε ε ε ∞ − ∆ = + ∆ La méthode ADE permet d’obtenir une relation entre D et E pour des modèles même non linéaires (ce que ne n’autorise pas la méthode RC à cause du problème lié à la définition CHAPITRE IV : Méthode des différences finies dans le domaine temporel 104 d’une relation de récurrence sur la sommation des champs électriques). Cependant, l’intérêt d’une telle méthode est essentiellement lié à l’aisance avec laquelle la relation entre D et E est exprimée dans le domaine temporel et par suite discrétisé en différence finie. IV.9 Conclusion Dans ce chapitre, nous avons présenté les fondements théoriques de la méthode des différences finies dans le domaine temporel « FDTD » et leur formulation numérique adaptée à l’analyse des structures à bande interdite photonique ainsi que leur implémentation informatique. Cette méthode de modélisation nous a permis de connaître la réponse spectrale et de calculer les distributions de champ dans des structures photoniques. Ensuite Nous avons exposé les différentes sources d’excitation dont l’impulsion est l’excitation la mieux adaptée. Enfin, nous avons présenté les conditions d’absorption aux limites et décrit la technique d’implémentation des PMLs, les conditions d’absorption les plus performantes aujourd’hui, dans l’algorithme de la méthode FDTD. Les matériaux dispersifs, représentés par le modèle de Debye et celui de Lorentz, ont été implémentés par la méthode de convolution récursive. CHAPITRE IV : Méthode des différences finies dans le domaine temporel 105 Bibliographie [1] K. S. Yee, Numerical solution of initial boundary value problems involving Maxwell's equations in isotropic media, IEEE Transactions on Antennas and Propagation, vol. 14, no. 3, pp. 302{307, 1966. [2] A.Taflove and M.E.Brodwin, IEEE Transactions on Microwave theory and Technique, MTT-23, No 8, August 1975. [3] A.Taflove and M.E.Brodwin, IEEE Transactions on Microwave theory and Technique, MTT-23, No 11, November 1975. [4] K. S. Kunz and K. M. Lee, IEEE Transactions on Electromagnetic Compatibility, EMC-20, No 2, May 1978. [5] A. Taflove and S. C. Hagness, Computational Electrodynamics, the Finite-Difference Time-Domain Method, 2nd ed. (Artech House, Norwood, MA, 2000). [6] R. Courant, K. Friedrichs et H. Lewy, Über die partiellen Differenzengleichungen der mathematischen Physik , Mathematische Annalen, vol. 100, no. 1, p. 32–74, 1928. [7] K. S. Kunz et R. J. Luebbers, The FDTD Method for Electromagnetics (CRC, Boca Raton, 1993). [8] K.R. Umashankar, et A.Taflove, IEEE Trans. Electromagnetic Compatibility, 24, 397(1982). [9] B. Engquist, A. Majda, Absorbing boundary conditions for the numerical simulation of waves, Mathematics of computation, vol.31,pp. 629-651, 1977. [10] A. Bayliss, E. Turkel, Radiation boundary conditions for wave-like equations, Comm.Pure Appl. Math.,vol. 23, pp. 707-725, 1980. [11] G.Mur, Absorbing boundary conditions for the finite difference domain electromagnetic field equations, IEEE Transaction Electromagnetic Compatibility, vol. 23, pp. 377-382, 1981. [12] R. Holland, L. Simpson, A free field EMP coupling and scattering code, IEEE, Transaction on nuclear sciences, vol. 6,n)24, pp. 2416-2421, 1977. [13] J.P. Berenger, A perfectly matched layer for the absorption of electromagnetic waves Journal of Computational physics, Vol. 114, pp. 185-200 (1994). [14] W.C. Chew, W.H.Weedon, A 3d perfectly matched medium from modified Maxwell’s equations with stretched coordinates, Microwave Opt. Technol. Lett 7, pp. 599-604, 1994. CHAPITRE IV : Méthode des différences finies dans le domaine temporel 106 [15] M. Kuzuoglu, R. Mittra, Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers, IEEE Microwave Guide Wave Lett 6, pp. 447-449, 1996. [16] J. Alan Roden, Stephen D. Gedney, Convolution PML (CPML) :an efficient FDTD implementation of the CFS-PML for arbitrary media, Microwave Opt Technol Lett, Vol 27, pp. 334-339, 2000. [17] D. Sullivan, Electromagnetic simulation using the FDTD method, (Wiley-IEEE Press, 2000). [18] J.P. Berenger, Perfectly Matched Layer (PML) for Computational Electromagnetics, (Morgan & Claypool, 2007). [19] Luebbers R.J., et al., 1990, A frequency-dependent finite-difference time-domain formulation for dispersive materials, IEEE Transactions on Electromagnetic Compatibility, 32, pp 222-227. [20] M. Okoniewskiand E. Okoniewska, Drude dispersion in ADE FDTD revisited, Electron. Lett., vol. 42, no. 9, pp. 503-504, 2006. CHAPITRE V Méthode de décomposition en ondes planes CHAPITRE V : Méthode de décomposition en ondes planes 108 V.1 Introduction La méthode de décomposition en ondes planes PWE (Plane Wave Expansion) s'est imposée comme l'un des outils de modélisation privilégié des cristaux photoniques [1, 2, 3] et figure, par ailleurs, parmi les premiers formalismes à avoir été employés afin de mettre théoriquement en évidence l'existence de bandes interdites pour les ondes électromagnétiques. Elle permet de représenter de façon assez directe, du point de vue du formalisme mathématique comme de la mise en œuvre numérique, la propagation de champs (électromagnétiques ou de déplacement) dans un milieu périodique. La technique PWE consiste à résoudre, dans l'espace fréquentiel, l'équation d'onde linéaire en développant le champ électromagnétique sur une base d'ondes planes. La méthode de décomposition en ondes planes est très efficace pour calculer les diagrammes de bandes de cristaux photoniques parfaitement périodiques. Elle permet de déterminer la fréquence, la polarisation, la symétrie et la distribution du champ pour les modes d'une structure photonique. Elle peut être adaptée pour étudier certaines structures non périodiques comme les guides d'onde ou les cavités grâce à la technique des supercellules [4] ; toutefois pour assurer une convergence des calculs, cette méthode nécessite un nombre important d'ondes planes, ce qui va induire un temps de calcul élevé et limiter cette méthode. V.2 Equation de Helmholtz Rappelons que tout phénomène électromagnétique est gouverné par les équations des Maxwell. Ces dernières conduisent à une équation d’onde (dite équation maîtresse ou encore équation de Helmholtz) qui, dans un milieu linéaire, isotrope, non magnétique (perméabilité magnétique relative égal à 1) de constante diélectrique complexe ε et en absence de sources, s’écrit (pour le champ magnétique) sous la forme: 2 1 ( ) ( ) ( ) H r H r r c ω ε | | | | ∇× ∇× = | | \ . \ . (5.1) CHAPITRE V : Méthode de décomposition en ondes planes 109 Par conséquent, les équations de Maxwell à travers l’équation de Helmholtz, sont représentées dans le domaine fréquentiel et transformées en un problème aux valeurs propres dont la résolution permet d’obtenir les relations de dispersion reliant la fréquence au vecteur d’onde [5]. Dans ce cas, ( ) H r est une fonction propre d’une certaine structure périodique définie par une fonction diélectrique périodique ( ) r ε et associée à la valeur propre ( ) 2 / c ω . La résolution analytique de l’équation (5.1) s'avère alors impossible. La méthode PWE est l'un des outils numériques destinés à résoudre cette équation. V.3 Structure de bandes des cristaux photoniques unidimensionnels V.3.1 Position du problème Le problème du calcul de structure de bandes d'un cristal photonique est l'obtention de la relation de dispersion qui représente la dépendance des fréquences propres ou les fréquences de résonance de cette structure avec le vecteur d’onde. Il s’agit de trouver des fonctions propres d'une structure périodique spécifique et des valeurs propres correspondantes et de déterminer ensuite la forme de la relation de dispersion. Afin d'obtenir cette relation, il est nécessaire de résoudre le problème aux valeurs propres formulé pour l'équation de Helmholtz, à l'intérieur de la structure périodique infinie. Pour un cristal photonique à une dimension, l’équation (5.1) devient [5]: 2 ˆ ( ) ( ) H x H x c ω | | Θ = | \ . , 1 ˆ ( ) x x x ε ∂ ∂ Θ = − ∂ ∂ (5.2) Etant donnée la périodicité de la constante diélectrique, on peut appliquer le théorème de Bloch à l’équation (5.2) pour développer le champ H en ondes planes, en exploitant l’analogie avec l'équation de propagation d'une onde électronique dans un cristal ordinaire où le potentiel périodique est donné par l'arrangement régulier des atomes. Le champ magnétique peut alors prendre la forme [6]: ( ) ( ) exp( ) x H x h x ik = (5.3) où h(x) est une fonction vectorielle périodique telle que : ( ) ( ) h x h x T = + . CHAPITRE V : Méthode de décomposition en ondes planes 110 En injectant (5.3) dans l’équation de Helmholtz (5.2), on arrive à l’expression suivante : 2 2 1 ( ) exp( ) ( ) exp( ) ( ) x x h x ik h x ik x x x c ω ε | | ∂ ∂ − ⋅ = ⋅ | ∂ ∂ \ . (5.4) La périodicité de la permittivité diélectrique ( ) x ε rend la solution d’un tel problème beaucoup plus complexe que le cas d’un milieu uniforme et, par conséquent, la relation de dispersion prendra une forme plus complexe. La forme typique de la relation de dispersion d’un cristal photonique unidimensionnel par rapport à la relation de dispersion d’un milieu uniforme est donnée dans la figure 5.1 Figure 5.1 : Diagramme de dispersion d’un milieu périodique. V.3.2 Calcul de structure de bandes Afin de résoudre l'équation aux valeurs propres (4.18), on devrait employer quelques méthodes approximatives qui utilisent la périodicité de la distribution de constante diélectrique. Par exemple, nous savons que la fonction propre d'une structure périodique infinie devrait également être infinie et périodique. C'est pourquoi le théorème de Bloch Milieu uniforme Milieu périodique CHAPITRE V : Méthode de décomposition en ondes planes 111 devrait être employé pour la représentation de la fonction propre du cristal photonique. Le champ H(x) satisfait le théorème de Bloch et peut être décomposé sous la forme [6]: , ( ) ( ) exp( ) k n H x h x jkx = ⋅ (5.5) où , ( ) k n h x est une fonction périodique de même périodicité que le réseau, k désigne le vecteur d'onde de Bloch et n le numéro de bande. Ainsi, en présence des fonctions infinies, il n'est pas possible de poursuivre le calcul. Il est donc commode de développer la fonction (5.5) en séries de Fourier : , ( ) ( ) exp( ( ) ) k n G H x h G j k G x = ⋅ + ∑ (5.6) De même pour l’inverse de la fonction diélectrique ( ) x ε : 1 ( ) exp( ) ( ) G G G jG x x χ ε ′′∈ ′′ ′′ = ⋅ ∑ (5.7) où k est le vecteur d’onde appartenant à la première zone de Brillouin, G le vecteur du réseau réciproque et , ( ) k n h G et ( ) G χ ′′ sont les coefficients de Fourier. Après le développement de toutes les fonctions infinies, nous les substituons dans l'équation de Helmholtz (5.2) : , 2 , 2 ( ) exp( ) ( ) exp( ( ) ) ( ) exp( ( ) ) 0 k n G G k n G G jG x h G j k G x x x h G j k G x c χ ω ′′ ′ ∂ ∂ ′′ ′′ ′ ′ ⋅ ⋅ + ∂ ∂ + ⋅ + = ∑ ∑ ∑ (5.8) et en tenant compte de l’égalitéG G G ′ ′′′ = + , l'expression ci-dessus devient [7, 8] : CHAPITRE V : Méthode de décomposition en ondes planes 112 , 2 , 2 ( ) exp( ( ) ) ( ) exp( ( ) ) ( ) exp( ( ) ) 0 k n G G k n G G G j G G x h G j k G x x x h G j k G x c χ ω ′ ∂ ∂ ′ ′ ′ ′ − ⋅ − ⋅ ⋅ ⋅ + ∂ ∂ + ⋅ + = ∑∑ ∑ (5.9) En dérivant et en combinant les exposants, on obtient : , 2 , 2 ( ) ( ) exp( ( ) ) [ ( ) ( )] ( ) exp( ( ) ) 0 k n G G k n G G G h G j k G x j k G j k G h G j k G x c χ ω ′ ′ ′ ′ − ⋅ ⋅ + × + ⋅ + + ⋅ + = ∑∑ ∑ (5.10) et la projection sur la base exp( ( ) ) j k G x + donne l’équation maîtresse d’un cristal photonique unidimensionnel : 2 , , 2 ( )[( ) ( )] ( ) ( ) 0 k n k n G G G k G k G h G h G c ω χ ′ ′ ′ ′ − − + ⋅ + ⋅ + = ∑ (5.11) L’équation (5.11) représente un système linéaire de dimension infinie car il existe une infinité de vecteurs du réseau réciproque et les informations sur la distribution de constante diélectrique à l'intérieur de la structure photonique sont maintenant fournies sous forme de coefficients de Fourier. L'opérateur différentiel dans l'équation (4.11) est présenté sous forme d’une matrice dont l'élément peut être déterminé à partir de l'expression suivante : , ( )(( ) ( )) G G G G k G k G θ χ ′ ′ ′ = − + ⋅ + (5.12) L’ensemble de solutions du système d'équations (4.11) est alors l’ensemble des valeurs propres de l'opérateur différentiel matriciel dont la forme est comme suit: CHAPITRE V : Méthode de décomposition en ondes planes 113 1 1 2 1 1 1 2 2 2 2 1 2 ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ ˆ N N N N N N G G G G G G G G G G G G G G G G G G θ θ θ θ θ θ θ θ θ θ ′ ′ ′ ′ ′ ′ ′ ′ ′ = (5.13) La particularité principale de la matrice (5.13) est qu'elle est hermitienne [9]. La matrice hermitienne est celle avec les éléments qui remplissent la condition suivante : * , , ˆ ˆ G G G G θ θ ′ ′ = (5.14) V.3.3 Solution du problème aux valeurs propres La solution du problème aux valeurs propres pour un cristal photonique unidimensionnel est généralement représentée sous forme de structure de bandes : en abscisses, le vecteur d’onde k ; en ordonnées la fréquence normalisée / 2 a c ω π . Dans le cas considéré, la fréquence est normalisée par la période du cristal photonique et ne dépend pas de l'échelle de la structure. Figure 5.2 : Structure de bandes typique d’un cristal photonique 1d Les fréquences propres du cristal photonique commencent à partir de la fréquence nulle à un point k = 0. Plus haut, sur l'axe de fréquences la bande interdite photonique apparaît. CHAPITRE V : Méthode de décomposition en ondes planes 114 Dans le cas d’un cristal photonique unidimensionnel, le PBG apparaît entre presque chaque paire des bandes. V.3.4 Algorithme de la méthode PWE Le processus entier de calcul effectué, pour un cristal photonique unidimensionnel, peut être décrit par les opérations suivantes : Figure 5.3 : Processus de calcul par la méthode PWE V.4 Structure de bandes des cristaux photoniques 2D et 3D V.4.1 Cas d'un cristal photonique 3D Comme dans le cas du cristal photonique unidimensionnel, le calcul de la structure de bandes exige l’équation de Helmholtz (5.1) obtenue à partir des équations de Maxwell. Dans une structure photonique 3D, la variation de la constante diélectrique est périodique dans toutes les directions possibles. Ecriture de l’expression pour le calcul des coefficients de Fourier. Limitation de la variation du vecteur d'onde dans la zone de Brillouin / / k T T π π = − Ensembles de G et de G’ varient entre les limites 2 / 2 / N T N T π π − où 2N + 1 est le nombre d'ondes planes considéré Ecriture de l'opérateur différentiel matriciel pour chaque vecteur d'onde dans l’intervalle choisie et calculer les états propres de la matrice obtenue. CHAPITRE V : Méthode de décomposition en ondes planes 115 Puisque nous recherchons des états propres d’une structure périodique infinie, la distribution spatiale des composants E et H du champ électromagnétique peut être représentée sous forme de fonctions de Bloch [6]: , ( ) ( ) jk r k n H r H r e ⋅ = ⋅ (5.15) où , ( ) k n H r est une fonction périodique de même périodicité que le réseau. Les indices indiquent que la fonction périodique est différente pour chaque vecteur d'onde k et nombre d'état propre n. la fonction , ( ) k n H r doit satisfaire la condition suivante : , , ( ) ( ) k n k n H r R H r + = (5.16) Comme dans le cas d’un cristal photonique unidimensionnel, la solution directe de l'équation de Helmholtz (5.1) pour une structure périodique infinie est impossible. Cependant, la fonction , ( ) k n H r étant périodique, elle peut être développée en séries de Fourier : , , ( ) ( ) exp( ( ) ) k n k n G H r H G j k G r ′ = + ⋅ ∑ (5.17) où G est un vecteur du réseau réciproque. Due à la périodicité de la fonction diélectrique ( ) r ε , il est également commode de développer l’inverse de cette fonction en séries de Fourier: 1 ( ) exp( ) ( ) G G jG r r χ ε = ⋅ ⋅ ∑ (5.18) où χ (G) sont les coefficients de Fourier. CHAPITRE V : Méthode de décomposition en ondes planes 116 En utilisant les fonctions de base exp( ) G r ⋅ , on peut transformer l’équation d’onde à la représentation de vecteurs d'onde. En remplaçant dans l'équation (5.1), les expressions (5.17) et (5.18), une équation aux valeurs propres est obtenue : { } 2 , 2 , , ( )( ) ( ) ( ) ( ) k n k n k n G G G k G k G H G H G c ω χ ′ ′ ′ ′ − − + × + × = ∑ (5.19) Cette équation représente « l’équation maîtresse » d’un cristal photonique 3D et ses solutions représentent les états propres de cette structure. L’équation (5.19) représente un système linéaire de dimension infinie car il y a une infinité de vecteurs G du réseau réciproque. La diagonalisation, qui doit être effectuée pour chaque valeur de k , permet alors de déterminer les valeurs propres , k n ω (n servant à numéroter les valeurs propres). Les valeurs de k sont limitées à certaines directions de symétrie de la première zone de Brillouin. Les courbes de dispersion du cristal photonique sont alors obtenues. Elles représentent les diagrammes de bandes du cristal. D’une manière générale, quand les vecteurs k décrivent la première zone de Brillouin, les fréquences , k n ω recouvrent continûment le spectre d’énergie. Cependant, dans certains cas, il existe des domaines d’énergie dans lesquels aucun mode , k n ω n’est accessible, ce sont les bandes interdites photoniques. V.4.1 Cas d'un cristal photonique bidimensionnel Pour passer à un problème 2D [10], il faut commencer par faire l'hypothèse que le problème est invariant par translation dans une direction. On peut alors décomposer les équations de Maxwell en deux systèmes indépendants : Transverse Electrique TE et Transverse Magnétique TM, faisant tous les deux intervenir l'équation de Helmholtz scalaire. Les équations de Helmholtz obtenues, pour les polarisations TE et TM à partir des équations de Maxwell, présentent des différences par rapport à ceux en 3D. CHAPITRE V : Méthode de décomposition en ondes planes 117 Les équations de Helmholtz pour le cas 2D prennent la forme suivante: 2 2 1 1 ( ) ( ) ( ) ( ) z z H r H r x r x y r y c ω ε ε ¦ ¹ ∂ ∂ ∂ ∂ ¦ ¦ − + = ´ ` ∂ ∂ ∂ ∂ ¦ ¦ ¹ ) (5.20) où r est un vecteur dans le plan (x,y) et la permittivité varie suivant z. Le processus de calcul de la structure de bandes d’un cristal bidimensionnel repose sur les mêmes étapes de calcul que dans le cas avec les structures photoniques 1D et 3D. Dans ce cas, l’équation aux valeurs propres pour les coefficients de Fourier prend la forme : 2 , 2 , , , , ( )( )( ) ( ) ( ) k n z k n z k n G G G k G k G H G H G c ω χ ′ ′ ′ − + + = ∑ (5.21) où G et G′ sont des vecteurs « in-plane » du réseau réciproque, k est le vecteur d’onde «in- plane » et , k n ω représente les fréquences propres des polarisations TE et TM (dans le cas d’un cristal photonique 2D, les valeurs de ces fréquences sont différentes). V.5 Développement de Fourier de la fonction diélectrique L’expression analytique pour les coefficients de développement de Fourier peut seulement être obtenue pour les cas les plus simples. En cas de présence des irrégularités dans le cristal photonique ou quand la forme des éléments de ce cristal diffère des formes standards (sphère, tige,…), il est nécessaire d'utiliser des méthodes numériques pour calculer les coefficients de Fourier. La méthode est basée sur la discrétisation de l'espace de la cellule unitaire et le remplacement de l'intégrale par une sommation [11]. La meilleure façon de discrétisation est de diviser la cellule unitaire par le maillage uniforme carré (un maillage avec des cellules identiques). Le processus de discrétisation est représenté sur la figure 5.4. CHAPITRE V : Méthode de décomposition en ondes planes 118 Figure 5.4 : Etapes de discrétisation d’une cellule unitaire par le maillage uniforme carré Le maillage doit être assez dense pour reproduire la forme des frontières des éléments constituant la structure photonique. V.6 Structure de bandes «off-plane » d’un cristal photonique 2D Dans le cas d’un cristal photonique bidimensionnel, la propagation hors plan « off- plane » est caractérisée par une constante de propagation z k non nulle suivant z. Donc, il n’est plus possible de séparer le système en deux sous-systèmes comme auparavant. Les cas TE et TM se trouvent donc mélangés et ne peuvent pas être traités séparément. Figure 5.5 : Propagation « off-plane » dans une structure photonique 2D Le calcul de la structure de bandes off-plane exige la modification de la zone de Brillouin. La modification consiste en une transition parallèle de la zone de Brillouin le long CHAPITRE V : Méthode de décomposition en ondes planes 119 de la direction z par la valeur de la composante du vecteur d’onde z k et mène à la modification d'opérateur différentiel matriciel : , , , , , , , , , ˆ ( ) ( ) ( ) x y x y x y z x y x y z x y G G H G G k G k G χ ′ ( ′ ′ = − + ⋅ + ¸ ¸ (5.22) où , , x y z x y z k k k k = + + est le vecteur d’onde en 3D. Figure 5.6 : Structure de bandes a) in-plane et b) off-plane V.7 Structure de bandes d’un cristal photonique avec défaut Lorsqu’un défaut est introduit dans la périodicité d’un cristal photonique possédant une bande interdite, de nouveaux modes permis pour le champ électromagnétique peuvent apparaître pour des fréquences se trouvant dans le gap photonique. Le calcul de la structure de bandes de ces objets permettrait de connaître la position et la dispersion de ces modes [9, 12]. Dès que des défauts sont introduits, la périodicité des structures est rompue et la méthode du développement en ondes planes ne peut plus s’appliquer. Pour continuer à utiliser cette méthode, une nouvelle périodicité peut être introduite artificiellement : le défaut est placé au centre d’une cellule de base comprenant plusieurs rangées du réseau d’origine que l’on appelle « supercellule ». Cette dernière est ensuite répétée suivant les directions du réseau de base pour donner naissance à un nouveau réseau parfaitement périodique. CHAPITRE V : Méthode de décomposition en ondes planes 120 Figure 5.7 : Exemple d’une nouvelle périodicité « supercellule » Cette méthode est appelée la méthode des supercellules. L’approximation que constitue cette méthode dépend directement du couplage entre les "défauts" du réseau constitué de supercellules. Il est évident que, plus les défauts sont éloignés, plus le couplage sera faible et plus les propriétés dispersives de la structure constituée de supercellules seront proches de celles de la structure ne comprenant qu’un défaut [13]. Figure 5.8 : Structure de bandes d’un cristal photonique a) sans défaut b) avec défaut La méthode de la supercellule s'applique aussi au calcul du diagramme de bandes de cristaux photoniques en présence de défauts étendus. CHAPITRE V : Méthode de décomposition en ondes planes 121 V.8 Conclusion La résolution analytique de l’équation maîtresse (chapitre 2), qui régit la propagation des ondes électromagnétiques dans les cristaux photoniques, n’est possible que dans des cas extrêmement simples. Pour les structures photoniques réelles, des méthodes de résolutions numériques sont nécessaires. Dans ce chapitre nous avons présenté la formulation et l’algorithme de base de la méthode de décomposition en ondes planes PWE qui sert à résoudre numériquement et de manière très efficace cette équation. Nous avons développé les fonctions périodiques, obtenues en appliquant le théorème de Bloch, en séries de Fourier pour pouvoir transformer la résolution des équations de Maxwell en un problème classique de diagonalisation de matrice. La résolution du problème aux valeurs propres obtenu permet de déterminer le spectre de fréquences propres et, par conséquent, de construire la structure de bandes du cristal photonique étudié. CHAPITRE V : Méthode de décomposition en ondes planes 122 Bibliographie [1] K. M. Ho, C. T. Chan, and C. M. Soukoulis, Existence of a photonic gap in periodic dielectric structures, Phys. Rev. Lett., vol. 65, no. 25, pp. 31523155, 1990. [2] M. Plihal and A. A. Maradudin, Photonic band structure of twodimensional systems: The triangular lattice, Phys. Rev. B, vol. 44, no. 16, pp. 85658571, 1991. [3] P. R. Villeneuve and M. Piché, Photoinc band gaps in twodimensional square and hexagonal lattices, Phys. Rev. B, vol. 46, no. 8, pp. 49694972, 1992. [4] R. D. Meade, K. D. Brommer, A. M. Rappe, and J.D. Joannopoulos, Existence of a photonic band gap in two dimensions, Appl. Phys. Lett., vol. 61, no. 4, pp. 495497, 1992. [5] J. D. Joannapolous, R. D. Meade, and J. N. Winn, Photonic Crystals - Molding the Flow of Light, Princeton University Press, Second edition (2008). [6] K. Sakoda, Optical Properties of Photonic Crystals (Second Edition). (Springer, Heidelberg, 2005). [7] J.M. Lourtioz, H. Benisty, V. Berger, et al., Photonic Crystals, Towards Nanoscale Photonic Devices. (Springer, Heidelberg, 2005). [8] K. Yitzhak, An introduction to Harmonic Analysis, 2nd corrected ed. (Dover, NY, 1976). [9] H. S. Sözüer and J. W. Haus, Photonic bands: Convergence problems with the planewave method, Phys. Rev. B, vol. 45, no. 24, pp. 1396213972, 1992. [10] Maksim Skorobogatiy, Jianke Yang, Fundamentals of Photonic Crystal Guiding, (Cambridge University Press, 2009). [11] J.Igor A. Sukhoivanov, Igor V. Guryev, Photonic crystals: Physics and Practical Modeling,. (Springer, Heidelberg, 2009). [12] S. G. Johnson and J. D. Joannopoulos, Photonic Crystals: The Road from Theory to Practice. Boston, MA: Kluwer Academic Publishers, 2002. [13] A. Yariv, Y. Xu, R. K. Lee, and A. Scherer Coupled-resonator optical waveguide: a proposal and analysis, Optics Lett 24, 711 (1999). CHAPITRE VI Conception et développement d'un logiciel de simulation CHAPITRE VI : Conception et développement d'un logiciel de simulation 124 VI.1 Introduction Un logiciel de simulation est un effort multidisciplinaire et la modélisation numérique des caractéristiques des cristaux photoniques nécessite en général un gros investissement en programmation et en analyse numérique et constitue une activité exigeante en compétence et en temps. Le développement d’un logiciel « maison » adapté aux besoins spécifiques de la recherche dans le domaine des matériaux à bande interdite photonique dont l’avantage principal est de pouvoir intégrer facilement et rapidement les dernières avancées, semble être un choix raisonnable. Notre objectif était la construction d’un environnement de simulation performant et riche en outils d’analyse, d’interfaces de représentation numérique et géométrique et de modules de visualisation. Cet environnement peut être aussi considéré comme une voie vers le développement d’un vrai « laboratoire virtuel » qui peut offrir aux physiciens et aux chimistes intéressés par le domaine des cristaux photoniques, l’opportunité de réaliser leurs expériences numériquement et de les assister pendant les procédés de synthèse et de configuration des matériaux à bande interdite photonique. VI.2 Description et architecture du logiciel Ce logiciel, nommé actuellement « PhcLab » (Photonic Crystals Laboratory) et destiné à l’analyse des matériaux à bande interdite photonique, est un environnement graphique de simulation basé sur les deux méthodes populaires : la méthode de différences finies dans le domaine temporel « FDTD » et celle de développement en ondes planes « PWE » détaillées précédemment. Il a été réalisé sous l’environnement de programmation C++Builder [1]. Nous l'avons choisi car il permet de conserver la rapidité d'exécution du C/C++ tout en simplifiant le processus de création de l'interface graphique. Afin de faciliter la saisie des paramètres de simulation, la représentation géométrique des structures photoniques et la visualisation des résultats, une interface graphique conviviale et flexible a été développée. La structure hiérarchique modulaire et l’exploitation des possibilités offertes par le concept orienté objet [2] assurent une plus grande souplesse de réutilisation des codes. Elles donnent au programme une meilleure extensibilité qui permet d'intégrer plus facilement les modifications et les améliorations ultérieures. CHAPITRE VI : Conception et développement d'un logiciel de simulation 125 La structure globale du logiciel est mise en évidence dans l’organigramme (fig. 6.1), qui décrit l’enchaînement des principaux modules ainsi que leurs plus importantes fonctions. Figure 6.1 : Organigramme global du logiciel. Module de simulation Simulateur FDTD • Traiter et déterminer les coefficients du calcul. • Mise à jour des champs. • Imposer les conditions aux limites absorbantes. Simulateur PWE • Calcul des coefficients de Fourier. • Calcul des modes propres. Module d’interface Windows • Routines pour créer l’application sous Windows. • Lire et analyser les données d'entrée. • Gérer l’interface graphique. Instructions et paramètres de simulation Module d’entrée • Propriétés des matériaux BIP. • Conditions aux limites absorbantes. • Sources d’excitation. • Paramètres du maillage. Module de sortie • Diagramme de dispersion • Structure de bandes. • Distribution du champ. • Réflectance et transmittance. Résultats de simulation CHAPITRE VI : Conception et développement d'un logiciel de simulation 126 VI.2.1 Module d’interface Windows Le module d’interface Windows se compose principalement des routines qui sont nécessaires pour créer l'application sous Windows et pour gérer l’interface graphique. Il permet de lire et analyser les données provenant du module d’entrée et de transmettre les instructions ainsi que les paramètres de simulation fournis par l’utilisateur, aux simulateurs. À la fin du calcul, ce module sert à la récupération des résultats de simulation. VI.2.2 Module d’entrée La fonction principale de ce module est de recevoir les informations fournis par l'utilisateur à travers l’interface graphique, la représentation géométrique de structures ou/et les fichiers de données. Ces informations représentent les paramètres de simulation, les paramètres de définition de structure et les instructions nécessaires pour contrôler la représentation graphique et numérique des résultats. VI.2.3 Module de sortie Ce module permet de représenter les résultats de simulation numériquement et graphiquement par l’intermédiaire des différents composants et outils de visualisation intégrés dans le logiciel. La distribution d’un champ électromagnétique dans une structure photonique et leur comportement dynamique, les valeurs de réflectance et de transmittance, le diagramme de dispersion et la structure de bandes sont généralement les résultats les plus importants qui sont pris en charge par le module de sortie. VI.2.4 Module de simulation Le module de simulation, qui représente le cœur du logiciel, comprend deux moteurs de calcul : le simulateur FDTD basé sur la méthode des différences finies dans le domaine temporel et le simulateur PWE basé sur la méthode de décomposition en ondes planes. VI.3 Simulateur FDTD Le simulateur FDTD dont l’organigramme est représenté sur la figure 6.2, se divise en trois parties ou étapes fonctionnelles: paramètres de structure et de simulation, génération du maillage et simulation FDTD. CHAPITRE VI : Conception et développement d'un logiciel de simulation 127 Figure 6.2 : Organigramme du simulateur FDTD. Au cours de la première étape, l'utilisateur spécifie les paramètres de simulation et les caractéristiques de la structure à analyser. Cette étape comprend les opérations suivantes : • La spécification de la gamme de fréquences, qui détermine la taille maximale autorisée de cellule et le nombre de cellules nécessaires pour simuler la structure ainsi que la largeur de bande du signal d'excitation. Étape de discrétisation spatiale Étape de discrétisation temporelle Étape de cartographie de matériaux Conception de la grille Démarrage du simulateur FDTD Initialisation des champs Étape de post-traitement Fin de simulation FDTD Résultats de simulation Mise à jour des composantes du champ H Mise à jour des sources d’excitation Mise à jour des composantes du champ E Application des conditions aux limites Incrémentation du temps par un pas temporel Dernière itération ? Maillage correct ? Non Oui Non Oui Définition de gamme de fréquences Définition de structure Définition de matériaux Spécification des sources d'excitation Spécification des conditions aux limites Définition des paramètres du maillage Informations valides ? Non Oui CHAPITRE VI : Conception et développement d'un logiciel de simulation 128 • La définition de la structure à analyser qui est établie soit en utilisant les outils géométriques disponibles dans le logiciel soit en choisissant un modèle prédéfini à partir de la bibliothèque. • La définition des matériaux qui peuvent être caractérisés à l’aide de trois paramètres : leur permittivité électrique, leur conductivité et leur perméabilité. • La spécification de la source utilisée pour exciter la structure photonique en indiquant généralement le type (gaussien, sinusoïdal, …), la largeur de bande, la fréquence et la position de cette source. • La sélection de types et de propriétés des conditions d’absorption aux limites choisies pour limiter le domaine du calcul. • La définition des paramètres du maillage utilisé par l’algorithme de génération du maillage uniforme. Après la configuration des paramètres de simulation et la définition de la structure photonique à analyser, le simulateur passe au niveau suivant qui consiste en la génération automatique du maillage. L'algorithme de génération du maillage effectue cette tâche en deux étapes : la première étape est la discrétisation spatiale et par conséquent la création d’une grille ; la deuxième est la formation d’une cartographie pour les matériaux utilisés, dont chaque cellule de la grille est assignée aux propriétés de matériau correspondant. Après avoir créé le maillage et déterminé la localisation et la taille de chaque cellule de la grille, la simulation FDTD peut être commencée. Avant que le procédé de progression temporelle (time-stepping procedure) ne commence, une étape d’initialisation des composantes des champs électriques et magnétiques et des coefficients de mise à jour de chaque cellule doit être accomplie. Les multiplicateurs constants et d'autres structures de données qui n'ont pas besoin d'être calculés ou assigné à chaque pas temporel sont également initialisés durant cette étape. Pendant l’étape de post-traitement, les quantités d'intérêt sont calculées. Après l'accomplissement de cette tâche, la simulation de FDTD est complète et les résultats sont rendus disponibles à l'utilisateur dans l'interface graphique. CHAPITRE VI : Conception et développement d'un logiciel de simulation 129 VI.4 Simulateur PWE Le simulateur PWE est résumé schématiquement par l’organigramme de la figure 6.3 et se compose de deux blocs fonctionnels: le bloc des paramètres de structure et de simulation et celui de simulation PWE. Figure 6.3 : Organigramme du simulateur PWE. Paramètres physiques de structure Limite la variation du vecteur d’onde k dans la zone de Brillouin Paramètres géométriques de structure Nombre d’ondes planes (2N+1) 2 Calcul des coefficients de Fourier Construire la matrice caractéristique Résoudre le problème aux valeurs propres Valeurs propres pour tout k dans la zone irréductible de Brillouin ? Groupe les valeurs propres dans des banes Informations valides ? Démarrage du simulateur PWE Fin de simulation PWE Résultats de simulation Étape de post-traitement Oui Non Oui Non CHAPITRE VI : Conception et développement d'un logiciel de simulation 130 Le premier bloc présente l’ensemble des paramètres requis pour lancer le calcul PWE. La structure photonique est définie par la période ou la constante du réseau, par les paramètres qui caractérisent la géométrie de la structure et les permittivités des matériaux. Le nombre d'ondes planes utilisées, qui doit être suffisamment grand (sans affecter sérieusement la convergence), détermine la précision du calcul. Avant de commencer le calcul de structure de bandes photoniques, les coefficients de Fourier sont calculés, soit analytiquement, soit numériquement à l’aide des algorithmes de Transformée de Fourier Rapide (Fast Fourier Transform, FFT). En spécifiant le vecteur d'onde k dans la zone de Brillouin irréductible qui représente la cellule élémentaire de l’espace réciproque, la matrice caractéristique est construite. Dans le cas d’un cristal photonique bidimensionnel, la matrice caractéristique de la polarisation TM est différente de celle de la polarisation TE. Le problème aux valeurs propres représenté par la matrice caractéristique est résolu par un algorithme standard de calculs matriciels. Ces dernières donnent une série de fréquences propres et de vecteurs propres correspondants. En répétant ce procédé, des modes propres sont trouvés pour un ensemble de vecteurs d'onde dans la zone de Brillouin irréductible. Ces modes sont énumérés de la plus basse à la plus haute fréquence. Une fois que nous avons la structure de bande dans la zone de Brillouin irréductible, le reste de la zone de Brillouin est obtenue par des opérations de symétrie. Après l'exécution de l’étape post-traitement, la simulation de PWE est complète et les résultats sont rendus disponibles à l'utilisateur dans l'interface graphique. L’ensemble des courbes de dispersion, qui forme le diagramme de bandes du cristal photonique, peuvent être affichées à l’aide d’un outil de visualisation. VI.5 Interface graphique L’interface graphique (Graphical User interface, GUI) du logiciel dont la fenêtre principale est montrée sur la figure 6.4 tente de satisfaire les critères de simplicité et de flexibilité. Le choix adopté porte sur une interface graphique principale définie par trois zones d’informations ou d’interactions. La première zone est la barre de menu, située en haut de la fenêtre principale, constitue l’élément majeur permettant de contrôler l’ensemble des fonctionnalités du logiciel. Cette barre permet de gérer les fichiers de données et les impressions. CHAPITRE VI : Conception et développement d'un logiciel de simulation 131 La deuxième zone consiste en une barre d’outils verticale dans laquelle sont placés des boutons qui assurent le déplacement rapide entre les catégories implémentés dans le logiciel. La liste des catégories comporte les structures unidimensionnelles, bidimensionnelles, tridimensionnelles et les structures de bandes. Figure 6.4 : Quelques fenêtres de l’interface graphique du logiciel. La zone principale de l’interface graphique se compose d’une boîte d’onglets proposant l’ensemble des fonctionnalités les plus couramment utilisées, tels que la saisie des paramètres, la représentation géométrique des structures, le paramétrage de l’espace de calcul CHAPITRE VI : Conception et développement d'un logiciel de simulation 132 et la visualisation des résultats. Cette boîte d’onglets permet une rapidité accrue quant à l’exécution des tâches qu’elle propose par rapport à un système de barre de menu. Les onglets disposent tous d’une barre de navigation rapide. VI.6 Validation du module de simulation L’étape de validation du logiciel de simulation est fondamentale. Cette opération consiste à effectuer des tests sur plusieurs structures photoniques et matériaux BIP. Ces tests seront portés sur le calcul des structures de bandes, par le simulateur PWE, de quelques cristaux photoniques et la simulation, par le simulateur FDTD, de la propagation de la lumière dans des structures photoniques unidimensionnelles et bidimensionnelles. VI.6.1 Validation du simulateur FDTD La figure 6.5 montre l’amplitude de Fourier d’une impulsion Gaussienne (500 Mhz) dans un matériau dispersif (le cas d’un cristal photonique unidimensionnel) décrit par le modèle de Debye et caractérisé par les paramètres suivants : ε =2, σ =0.01, 2 χ = , 0 0.001 t s µ = . Figure 6.5 : Amplitude de Fourier dans un milieu dispersif (Milieu de Debye). Les résultats de cette simulation (la distribution de l’amplitude, la valeur du coefficient de transmittance et la valeur de coefficient de réflectance) sont les mêmes que ceux rapportés par Sullivan [3]. L’introduction d’un défaut linéaire dans un cristal photonique conduit à l’apparition des modes localisés associés à ce défaut qui permettent le guidage des ondes électromagnétiques (fig. 6.6). CHAPITRE VI : Conception et développement d'un logiciel de simulation 133 Figure 6.6 : Propagation de la lumière dans un cristal photonique bidimensionnel en silicium (ε =11.6) avec un défaut linéaire. Notre calcul FDTD présente un très bon accord avec celui obtenu par Meep (MIT electromagnetic equation propagation) [4] ainsi que celui rapporté par Ozbay et al. [5]. Ces résultats comparatifs valident ainsi notre simulateur FDTD. VI.6.2 Validation du simulateur PWE La figure 6.7 représente la structure de bandes d’un cristal photonique en mode TE. Celui-ci est réalisé avec un réseau triangulaire sur l’Arséniure de Galium « GaAs » dont les paramètres sont : 13 ε = , le rapport r/a=0.48. Figure 6.7 : Structure de bandes d’un cristal photonique bidimensionnel en GaAs (polarisation TE, réseau triangulaire) Le diagramme de bandes montre l’existence d’une bande interdite totale dans l’intervalle 0.3600 - 0.5212, entre la première et la deuxième bande. Ce résultat est en accord avec celui rapporté par S. Guo [6]. CHAPITRE VI : Conception et développement d'un logiciel de simulation 134 La figure 6.8 montre le diagramme de bandes photoniques calculé par le simulateur PWE. La structure photonique est composée de trous d’air dans un milieu diélectrique, niobate de lithium d’indice n = 2:1421 avec un taux de remplissage de 0.2267%. Figure 6.8 : Structure de bandes d’un cristal photonique bidimensionnel en niobate de lithium (LiNbO3) (polarisation TE, réseau triangulaire) Dans ce cas (cas d’une structure triangulaire en polarisation TE), on note l’apparition d’une bande interdite photonique pour / 2 a c ω π compris entre 0.32 et 0.35. Ce résultat est en bon accord avec celui obtenu par MPB (MIT Photonic Bands) [7]. La figure 6.9 présente la structure de bandes d’un cristal photonique tridimensionnel en dioxyde de titane (TiO2) « Structure de type diamant ». Celle-ci correspond à un réseau CFC avec une constante diélectrique relativement faible ( ε =2.5). Figure 6.9 : Structure de bandes d’un cristal photonique tridimensionnel en TiO2 CHAPITRE VI : Conception et développement d'un logiciel de simulation 135 Le diagramme ne montre aucune bande interdite photonique à cause du faible indice de réfraction. Ce résultat est en accord avec celui obtenu par MPB (MIT Photonic Bands). VI.7 Conclusion Dans ce chapitre, nous avons tout d’abord présenté d’une manière générale les organigrammes et les algorithmes de base ainsi que les méthodes et les outils informatiques utilisées pour le développement du logiciel « PhcLab ». Ce logiciel à été développé en langage orienté objet C++Builder avec une structure modulaire et extensible et doté d’une interface graphique (GUI) conviviale. Les modules principaux du logiciel, notamment le simulateur FDTD basé sur la méthode de différences finies dans le domaine temporel et le simulateur PWE basé sur la méthode des ondes planes, ont été décrits. Par la suite, nous avons réalisé des tests de validation pour confirmer la fiabilité de chaque simulateur. Le simulateur FDTD a été validé par la simulation d’un milieu dispersif (milieu de Debye) et la propagation de la lumière dans un cristal photonique bidimensionnel en silicium avec un défaut linéaire. Le simulateur PWE a été validé par le calcul des structures de bandes des cristaux photoniques bidimensionnels en GaAs et LiNbO3 et une structure photonique de type diamant en TiO2. Les tests de validation ont donné des résultats en accord avec ceux rapportés dans la littérature et obtenus par Meep (MIT Electromagnetic Equation Propagation) et MPB (MIT Photonic Bands). CHAPITRE VI : Conception et développement d'un logiciel de simulation 136 Bibliographie [1] S.Kolachina, C++Builder 6 Developer’s Guide, (Wordware Publishing, 2002). [2] B.Bruegge and A. H. Dutoit, Object-Oriented Software Engineering: Conquering Complex and Changing Systems, (Prentice Hall, 1999). [3] D. Sullivan, Electromagnetic simulation using the FDTD method, (Wiley-IEEE Press, 2000). [4] Meep software is copyright of Massachusetts Institute of Technology (http://abinitio. mit.edu/wiki/index.php/Meep). [5] E. OZBAY, M. BAYINDIR, I. BULU, E. CUBUKCU, Investigation of localized coupled-cavity modes in two-dimensional Photonic Bandgap Structures, IEEE Journal of quantum electronics, vol. 38, N0. 7, July 2002. [6] S. Guo and S. Albin. Simple plane wave implementation for photonic crystal calculations. Optics express, 11(2) :167–175, 2003. [7] Johnson, S. and J. Joannopoulos, Block-iterative frequency-domain methods for Maxwell's equations in a planewave basis. Opt. Express, 2001. 8(3): p. 173-190. . CHAPITRE VII Applications des cristaux photoniques CHAPITRE VI I: Applications des cristaux photoniques 138 VII.1 Introduction Les propriétés de bande interdite, de diffraction, de réfraction, de localisation, de filtrage et de guidage des cristaux photoniques ne sont pas spécifiques d’un domaine de longueur d’onde particulier. Seuls les aspects de réalisation technique et les propriétés intrinsèques des matériaux qui composent le cristal changent évidemment, d’un domaine à l’autre, les échelles de longueur mises en jeu. Cette universalité explique la dynamique importante de recherche qui s’est créée autour des cristaux photoniques en rassemblant les physiciens, les opticiens, les chimistes et électroniciens. On retrouve ainsi la même diversité dans le champ des applications des cristaux photoniques avec une large couverture spectrale allant des ondes hyperfréquences jusqu’au domaine du visible. Les applications potentielles des cristaux photoniques sont très vastes et couvrent plusieurs domaines disciplinaires : réalisation des cavités résonantes de taille très réduite, des guides d’ondes [1], des lasers sans seuil [2], des filtres sélectifs [3], des multiplexeurs [4], des fibres optiques [5], des nouveaux composants optoélectroniques plus performant et compacts et des dispositifs reproduisant les principes opérationnels des différents composants d’un circuit intégré, en utilisant les photons comme porteur d’information à la place des électrons[6]. Ils trouvent également leurs applications dans le domaine de l'imagerie médicale [7, 8], cellules solaires à haut rendement [9], stockage d’information et d’énergie, développement de micro-capteurs chimiques et biologiques [10, 11, 12], blindage électromagnétique [13] et spectroscopie [14]. Parmi les applications potentielles des cristaux photoniques, nous ne présenterons dans ce chapitre que deux des applications d’intérêt chimique et biochimique : capteurs à base de cristaux photoniques et spectroscopie femtoseconde. CHAPITRE VI I: Applications des cristaux photoniques 139 VII.2 Capteurs à base de cristaux photoniques Les capteurs chimiques et biochimiques ont connu un développement croissant ces dernières années en raison de leur faible coût, de leur portabilité et de leurs nombreux domaines d’applications qui apparaissent aussi bien dans l’industrie automobile (contrôle des émissions de gaz), que l’industrie agroalimentaire (contrôle des procédés de fabrication), l’environnement ( détection des gaz toxiques), ou le biomédical. Parmi les différentes techniques de détection chimique, la détection optique offre un attrait de mesure en temps réel avec une très bonne sensibilité et une compatibilité avec les fibres optiques. L’application des cristaux photoniques en tant que capteurs constitue un domaine de recherche qui semble être très prometteur en raison de leur extrême miniaturisation, de leur haute sensibilité spectrale et de la possibilité de les intégrer aux MEMS [15]. Récemment, il y a eu plusieurs travaux de recherches utilisant les cristaux photoniques en tant qu’élément de détection en raison de leur structure de bande et du confinement de la lumière [16]. En exploitant la dépendance des propriétés optiques du cristal photonique aux caractéristiques physiques et géométriques du cristal lui-même, on peut distinguer les types des capteurs suivants : • Capteurs à base de cristaux photoniques avec défauts. • Capteurs à base de cristaux photoniques avec résonances guidées. • Capteurs à base d’opales artificielles. • Capteurs à cristaux photoniques basées sur les propriétés dispersives. • Capteurs à base de cristaux photoniques à fibres optiques. L’un des principes de fonctionnement de ces capteurs consiste en la mesure d’un changement d’indice de réfraction d’un élément sensible en fonction de la présence d’un analyte (substance à détecter). Généralement, la structure de tels capteurs se présente sous la forme d’un arrangement périodique de trous fonctionnalisés sur Niobate de Lithium (LiNbO3) ou Dioxyde de Titane (TiO2). Le choix de ces matériaux est motivé par les forts coefficients électro-optiques, acoustiques et non linéaires qui autorisent une commande rapide de l’indice de réfraction, permettant ainsi l’amélioration des performances du capteur. CHAPITRE VI I: Applications des cristaux photoniques 140 VII.2.1 Description de la structure La détection est effectuée par une monocouche sensible qui réagit avec l’analyte. La zone de détection est composée d’une structure photonique fonctionnant avec une monocouche sensible (fig. 7.1). Cette monocouche réagit avec l’analyte qui induit une variation de l’indice de réfraction, de l’épaisseur de la couche sensible et de son absorption. Cette structure doit être conçue pour être fortement sensible à la variation de l’indice de réfraction. Cette méthode peut être appliquée à la détection d’un grand nombre d’éléments chimiques et biologiques. Figure 7.1 : Exemple d’un capteur chimique à base d’un cristal photonique bidimensionnel en LiNbO3. La figure 7.1 présente un capteur chimique et biochimique [17] d’une structure à maille carrée. Les paramètres de ce capteur (notamment le rayon des trous et la période de la structure photonique) peuvent être déterminés en étudiant le diagramme de bandes du niobate de lithium de façon à bénéficier d’une structure simple (matrice sans défaut), avec la plus faible vitesse de groupe possible. VII.2.2 Modélisation de la structure L’étude de cette structure peut être réalisé a l’aide de la méthode de décomposions en ondes planes PWE (chap. 5) qui doit être modifié pour tenir compte des monocouches Monocouche sensible LiNbO3 CHAPITRE VI I: Applications des cristaux photoniques 141 sensibles [18]. L’étude consiste à calculer le diagramme de bandes et d’évaluer la variation de la structure de bandes en variant l’indice de réfraction. En utilisant la méthode FDTD (chap. 4), l’étude de cette structure consiste à calculer le spectre de transmission qui devrait donner un déplacement en longueur d’onde, d’une valeur déterminée, en présence de l’analyte. VII.3 Spectroscopie ultrarapide à base de cristaux photoniques La spectroscopie femtoseconde, qui utilise des sources émettant des impulsions ultracourtes, constitue un outil majeur en pleine évolution. Le degré de sophistication des expériences ne cesse de croitre y compris par le biais d’une ingénierie en amplitude et phase des impulsions ultracourtes. VII.3.1 Lasers à impulsions femtoseconde Les lasers à impulsions femtoseconde (1fs = 10 -15 Ces lasers sont également un outil permettant la manipulation à l’échelle moléculaire : c’est le domaine du « contrôle cohérent », dans lequel les impulsions peuvent être utilisées pour préparer le système afin d’en stimuler l’évolution dans une direction choisie. s) ont récemment ouvert de nouvelles perspectives d’application en analyse chimique, biochimique, traitement des matériaux, et physique de l’atmosphère. Leur avantage déterminant est de délivrer des puissances crêtes extrêmement élevées tout en ne déposant qu’une très faible énergie dans les échantillons. La disponibilité´ de sources laser d’impulsions toujours plus brèves permet désormais la résolution temporelle, qui donne accès au suivi détaillé de processus dynamiques inaccessibles jusqu’ici, tels que la dynamique interne, la relaxation liée aux interactions entre molécules ou atomes, les collisions réactives, et qu’il s’agit d’analyser et de comprendre. La production d'impulsion femtoseconde par une cavité laser se fait par la technique du blocage de mode (mode-locking). Le blocage de mode peut être actif (réalisé à l'aide d'un cristal électro-optique ou acousto-optique intra-cavité) ou passif. En pratique, la grande majorité des lasers femtoseconde utilisent la méthode du blocage de modes passif, qui permet d'obtenir des impulsions plus courtes. CHAPITRE VI I: Applications des cristaux photoniques 142 De manière générale, en optique guidée, il est possible de modifier les propriétés de dispersions par différentes techniques, en particulier celle des matériaux à cristaux photoniques. Plusieurs technologies de lasers femtoseconde utilisent ces possibilités. Les lasers basés sur la technologie des fibres à cristaux photoniques ont désormais atteint voire surpassé les performances offertes par les lasers solides conventionnels. Cette rupture a été possible grâce à l’exploitation de nouveaux régimes de propagation d’impulsions ultracourtes dans les fibres optiques à dispersion normale. Récemment, des impulsions femtoseconde de plus de 30 nJ d’énergie ont été générées dans un laser à fibre à bandes interdites photoniques 1D (fibre de Bragg) à cœur dopé à l’ytterbium [19]. VII.3.2 Description de la structure Une technique pour construit un spectromètre femtoseconde très compact [20] est basé sur l’utilisation d’un continuum de « lumière blanche » produit dans une fibre à cristaux photoniques. Ces fibres peuvent être conçues afin d’optimiser la largeur spectrale (450 – 1100 nm) du continuum tout en nécessitant des énergies modestes (nJ) telles que les oscillateurs Titane-Sapphire les produisent sans amplification coûteuse. La source laser du système est un oscillateur avec une cavité étendue, produisant des pulses de 40 fs pulses et de 15-18 nJ. La détection large bande se fait par une caméra CCD (Charge Coupled Device) qui peut acquérir jusqu’à 800 spectres/s, garantissant ainsi un excellent rapport signal à bruit. Figure 7.2 : section transverse de la fibre à bande interdite photonique. Une autre technique consiste en la génération d’impulsions femtoseconde dans un laser à fibre à bandes interdites photoniques tout solide à cœur dopé à l’ytterbium. Le verrouillage de modes est obtenu par la technique de l’évolution non-linéaire de polarisation assistée par un filtre spectral passif. Le laser émet des impulsions de 4 ps de durée avec un spectre de 12.5 nm de largeur centré autour de 1064 nm. Le laser délivre une puissance moyenne de 350mW à la cadence de 16.5 MHz, ce qui correspond à plus de 21 nJ d’énergie. Ces impulsions sont comprimées à l’extérieur de la cavité à 230 fs de durée. CHAPITRE VI I: Applications des cristaux photoniques 143 Bibliographie [1] R. D. Meade, A.Devenyi, J. D. Joannopoulos, O. L. Alerhand, D. A. Smith et K. Kash, Novel applications of photonic band gap materials: Low loss bends and Q cavities, Journal of applied physics 75, p. 4753(1994). [2] P. Pottier, C. Seassal, S. Letertre, J.L. Leclercq, P. Vik Torowitch, D. Cassagne, C. Jouanin, Triangular rand hexagonal high Q factor 2D photonic bandgap cavities on III- V suspended membranes, IEEE, Lightwave Technology, Journal of, November 1999, pp. 2058-2062 [3] M. Qiu, M. Mulot, M. Swillo, S. Anand, B. Jaskorzynak, A. Karlsson, M. Kamp, and A. Forchel, Photonic crystal optical filter based on contra-directional waveguide coupling, Applied Physics Letters Vol. 83, pp. 5121–52123 (2003). [4] Y. Akahane, M. Mochizuki, T. Asano, Y. Tanaka, and S. Noda, “Design of a channel drop filter by using a donor-type cavity with high-quality factor in a two-dimensional photonic crystal slab” Applied Physics Letters, Vol. 82, pp.1341–1343 (2003). [5] Temulkuran, B.; Hart, S. D.; Benoit, G.; Joannopoulos, J. D.; Fink, Y. Wavelength- Scalable Hollow Optical Fibers with Large Photonic Bandgaps for CO2 Laser Transmission Nature 2002, 420, 650. [6] X.J.M. Leijtens et M.K. Smit, S-matrix oriented CAD-tool for Photonic Integrated Circuits, Integrated Optic Devices II, Proceeding of SPIE, San Jose (USA), Janvier 1998. [7] L. Jylhä, I. Kolmakov, S. Maslovski, and S. Tretyakov, Journal of Applied Physics 99, 043102 (2006). [8] J. A. Schuller, R. Zia, T. Taubner, and M. L. Brongersma, Physical Review Letters 99, 107401 (2007). [9] K. C. Huang, M. L. Povinelli, and J. D. Joannopoulos, Applied Physics Letters 85, 543 (2004). [10] M. Loncar, A. Scherer, and Y. M. Qiu, Photonic crystal laser sources for chemical detection, Appl. Phys. Lett. 82, 4648-4650 (2003). [11] E. Chow, A. Grot, L. W. Mirkarimi, M. Sigalas, and G. Girolami, Ultracompact biochemical sensor built with two-dimensional photonic crystal microcavity, Opt. Lett. 29, 1093-1095 (2004). [12] F. Benabid, J. C. Knight, G. Antonopoulos et P. St. J. Russell, Stimulated Raman Scattering in Hydrogen-Filled Hollow-Core Photonic Crystal Fiber, Science, 298, 399 (2002). CHAPITRE VI I: Applications des cristaux photoniques 144 [13] J. Wang, J. Liang, H. Wu, W. Yuan, Y. Wen, Y. Song, L. Jiang, A facile method of shielding from UV damage by polymer photonic crystals, Poly. Inter. 0959-8103 (2008). [11] J. A. Schuller, R. Zia, T. Taubner, and M. L. Brongersma, Physical Review Letters 99, 107401 (2007). [15] Wonjoo Suh, M. F. Yanik, Olav Solgaard, and Shanhui Fan, Applied Physics Letters, 2003, Vol. 82, No 13, pp 1999-2001(2003). [16] T. Stomeo, M. Grande, A. Qualtieri, A.Passaseo, A.Salhi, M. Vittorio. Microelectronic Engineering, 2007, Vol. 84, issue 5-8, pp 1450-1453 (2007). [17] J. Dahdah, N. Courjal, F. Baïda, Structure photonique sur niobate de lithium dopé Mgo pour détection des gaz, Recueil des communications, pp. 213-215, Lannion, 20-22 octobre 2008. [18] K. Sakoda, Optical properties of photonic crystals, Springer Series in Optical Sciences, Vol. 80, Springer-Verlag, Berlin, (2005). [19] C. Lecaplain, A. Hideur, S. Février, P. Roy, Mode-locked Yb-doped Bragg fiber laser, Optics Letters, Vol. 34, no. 18, pp.2879-2881 (2009). [20] J. Léonard, N. Lecong, O. Crégut, S. Haacke, P. Vial, Ph. Leproux, V. Couderc, Broadband ultrafast spectroscopy using a photonic crystal fiber : application to the photophysics of malachite green, in preparation (2007). . Conclusion Générale et Perspectives Conclusion Générale et Perspectives 146 Dans le cadre de ce travail, nous nous sommes intéressés essentiellement à l’étude théorique des cristaux photoniques et le développement d’un environnement informatique de simulation, basé sur la méthode de différences finis dans le domaine temporel et la méthode de décomposition en ondes planes, pour analyser ces structures ainsi que des applications d’intérêt chimique et biochimique. Dans un premier temps, nous nous sommes focalisés sur l’étude théorique des matériaux à bande interdite photonique en exploitant l’analogie formelle qui existe entre les équations de Maxwell régissant la propagation des ondes électromagnétiques dans un milieu diélectrique et l’équation de Schrödinger vient de la périodicité géométrique du cristal atomique. A partir des formes fréquentielles des équations de Maxwell et en utilisant les outils développés en physique de l’état solide et en mécanique quantique, nous avons reproduit la forme générale (pour des structures en trois dimensions) de l’équation de Helmholtz. Nous avons montré que cette équation vectorielle représente un problème aux valeurs propres avec un operateur Hermitien et linéaire (pour le champ magnétique ( ) H r ) et admet par conséquent comme solutions des valeurs propres réelles. La loi d’échelle qui implique que le même comportement physique sera observé si l’on change la longueur d’onde et les dimensions du cristal dans les mêmes proportions à été montré. Ensuite, nous avons présenté les fondements théoriques de la méthode des différences finis dans le domaine temporel « FDTD » basé sur la discrétisation des équations de Maxwell et utilisée généralement pour la modélisation de la propagation des ondes électromagnétiques dans des milieux matériels. Nous avons formalisé cette méthode et implémenté des algorithmes pour modéliser les cristaux photoniques et connaître le comportement de la lumière (ou généralement les ondes électromagnétiques) dans ces structures périodiques. Les conditions d’absorption aux limites « PML » pour limiter l’espace du calcul et par conséquent la consommation des ressources informatiques ainsi que les modèles de Debye et de Lorentz pour les matériaux dispersifs ont été implémentés dans ces algorithmes. Pour résoudre l’équation de propagation obtenue à partir des équations de Maxwell représentées dans le domaine fréquentiel et construire ainsi le diagramme de dispersion, nous avons formulé et développé des algorithmes basés sur la méthode de décomposition en ondes planes « PWE ». Nous avons employé le théorème de Bloch et les séries de Fourier pour trouver la matrice caractéristique d’une structure photonique infinie dont les valeurs propres Conclusion Générale et Perspectives 147 sont les solutions de l’équation de propagation. Pour traiter le cas de présence des défauts intentionnels dans les structures photoniques, la technique de supercellules à été implémentée dans les codes PWE. Pour atteindre notre objectif, un logiciel de simulation et d’analyse des cristaux photoniques avec une interface graphique (GUI) et une structure modulaire a été réalisé. Ce logiciel, nommé actuellement « PhcLab » (Photonic Crystals Laboratory) et développé en langage orienté objet C++Builder, se compose de plusieurs modules dont le module principal comporte deux simulateurs reposent sur les méthodes FDTD et PWE. Plusieurs tests ont été effectués pour valider les codes développés et la performance des simulateurs. Nous avons utilisé le simulateur FDTD pour réaliser la simulation de la propagation d’une impulsion Gaussienne dans un cristal photonique unidimensionnel en matériau dispersif (modèle de Debye) et la propagation de la lumière dans un cristal bidimensionnel en silicium avec un défaut linéaire (guide d’onde). En suite nous avons calculé les structures de bandes à travers le simulateur PWE pour des cristaux photoniques bidimensionnels en GaAs et LiNbO3 et une structure photonique de type diamant en TiO2. Ces tests de validation ont donnés des résultats en accord avec ceux rapportés dans la littérature et obtenues par Meep (MIT Electromagnetic Equation Propagation) et MPB (MIT Photonic Bands). Le dernier volet a été réservé à la description de deux applications, à base de cristaux photoniques, d’intérêt chimique et biochimique. L’une d’entre elles concerne la réalisation des capteurs chimiques et biochimique à base de cristaux photoniques et l’autre concerne la spectroscopie femtoseconde. Dans le prolongement de ce travail, nous envisageons les perspectives suivantes : Etude théorique et modélisation des cristaux photoniques métallo-diélectriques, organiques et biologiques dont les constantes caractéristiques et les formes géométriques sont plus complexes ainsi que les métamatériaux dont l’indice de réfraction est négatif. Etude et modélisation des matériaux nanophotoniques et des cristaux photoniques magnéto-optiques ainsi que l’optique non linéaire dans les cristaux photoniques. Conclusion Générale et Perspectives 148 Amélioration des codes FDTD et PWE pour prendre en charge des nouvelles propriétés et des structures photoniques plus complexes et pour fournir des résultats plus précis sur une large bande de fréquences. Augmentation de la vitesse du calcul par optimisation et parallélisation des codes FDTD et par l’orientation des calculs vers le processeur graphique « GPU » qui est au moins six fois plus rapide que le microprocesseur. Améliorations au niveau de l’interface graphique du logiciel « PhcLab » tels que l’intégration d’un module CAD (computer Aided Design) plus efficace pour les structures photoniques 3D et d’une bibliothèque des structures photoniques personnalisables préconfigurés (notamment en 3D) ainsi que d’autre outils de visualisation. Développement de nouvelles applications d’intérêt chimique et biochimique à base de cristaux photoniques. Annexes Annexe I : Transformée de Fourier 150 I.1 Transformation de Fourier pour les fonctions intégrables La transformée de Fourier F est une opération qui transforme une fonction intégrable sur en une autre fonction, décrivant le spectre fréquentiel de cette dernière. Si f est une fonction intégrable sur , sa transformée de Fourier est la fonction ˆ ( ) F f f = donnée par : 1 ˆ ( ) : ( ) ( ) 2 i t F f f f t e dt ω ω ω π +∞ − −∞ = ∫ (1.1) En physique, la transformation de Fourier permet de déterminer le spectre d'un signal. On peut généraliser la définition de la transformée de Fourier à plusieurs variables, et même sur d'autres groupes que le groupe additif . I.2 Transformée de Fourier inverse Si la transformée de Fourier de f est elle-même une fonction intégrable, la formule dite de transformation de Fourier inverse, opération notée 1 F − , est celle qui permet (sous conditions appropriées) de retrouver f à partir des données fréquentielles : 1 ˆ ( ) ( ) 2 i t f t f e d ω ω ω π +∞ −∞ = ∫ (1.2) I.3 Transformée de Fourier rapide La transformée de Fourier rapide FFT (Fast Fourier Transform) est un algorithme de calcul de la transformée de Fourier discrète. Cet algorithme est couramment utilisé en traitement numérique du signal pour transformer des données discrètes du domaine temporel dans le domaine fréquentiel. Sa complexité varie en ( ln ) n n Ο avec le nombre de points n, alors que la complexité du calcul de base s'exprime en 2 ( ) n Ο . Ainsi, pour n=1024, le temps de calcul de l'algorithme rapide peut être 100 fois plus petit que le calcul utilisant la formule de définition de la transformée de Fourier discrète. Annexe I : Transformée de Fourier 151 Soient x 0 , ...., x n-1 2 1 0 i n jk n i k k f x e π − − = = ∑ des nombres complexes. La transformée de Fourier discrète est définie par la formule suivante : j=0,….,n-1 (1.3) ou en notation matricielle : 2 0 0 2 1 1 1 2 4 2( 1) 2 2 1 2( 1) ( 1) 1 1 1 1 1 1 1 1 1 n n n n n n n f x f x f x f x ω ω ω ω ω ω ω ω ω − − − − − − − = (1.4) Evaluer ces sommes directement coûte (n − 1) 2 produits complexes et n(n − 1) sommes complexes alors que seuls (n/2)(log 2 (n) − 2) produits et nlog 2 (n) sommes sont nécessaires avec la version rapide. Comme la transformée de Fourier inverse discrète est équivalente à la transformée de Fourier discrète, à un signe et facteur 1/n près, il est possible de générer la transformation inverse de la même manière pour la version rapide. Annexe II : Origine de la bande interdite photonique 152 Théorème variationnelle électromagnétique Le théorème variationnelle électromagnétique est l’analogue du principe variationnelle de la mécanique quantique. Afin de mieux comprendre le concept et l’origine de la bande interdite photonique, nous exploiterons deux propriétés générales des problèmes aux valeurs propres hermitiens. La première concerne l’orthogonalité des fonctions propres donnée par la relation suivante : ( )* ( ) 0 m n k k H H ⋅ = ∫ , avec m n ≠ (2.1) et la deuxième implique que la plus basse bande résout le problème variationnel : 2 2 1 2 ( ) / ( ) min k k ik H k c H ε ω ∇+ × = ∫ ∫ (2.2) La deuxième bande doit satisfaire le même problème variationnel, avec la contrainte supplémentaire qu’il soit orthogonale à la première bande selon (2.1), et ainsi de suite. Pour minimiser le problème variationnel (2.2), le champ électrique k E pour la plus basse bande devrait essayer de satisfaire deux conditions: leur rotationnel ne devrait pas être trop grande (il devrait varier lentement) et il devrait être concentré dans les régions où la constante diélectrique ε est haute. De même, la dixième bande veut également osciller lentement et être concentrée dans les zones à hauteε . Cependant, la formule (2.1) exige qu’elle doive être orthogonale à la première bande. Afin d'accomplir ceci, la deuxième bande doit être concentrée dans les zones à basse ε et/ou avoir une oscillation de signe opposé à l'intérieur des zones à haute ε pour avoir un intégral nul. Les deux possibilités encourront une grande augmentation de la fréquence dans l’équation (2.2). Il y aura généralement une grande différence entre les premières et les deuxièmes bandes et on obtient donc un « gap » à chaque point k . Une bande interdite omnidirectionnelle peut être obtenue par le recouvrement des gaps à tous les points k dans une certaine gamme de fréquences. Pour cette raison, un contraste assez grand de ε entre les deux régions diélectriques est généralement nécessaire pour rendre les gaps dans les différentes directions assez grandes pour assurer l’intersection.
Comments
Report "Etude théorique et simulation des cristaux photoniques"