Applied Science, Faculty of
Engineering, School of (Okanagan)
DSpace
UBCO
Warner, Cael
Master of Applied Science - MASc
University of British Columbia
Multiple electrodynamic theories have existed for over a century, which has led to controversy over a field-based description of light-matter interactions. Experiments that aim to support or disprove an electrodynamic theory benefit from accurate numerical simulation of each theory's predictions of light-matter interactions. However, most experiments and their numerical simulations are limited to cases of rigid-body translation or elastic deformation, which are equivalently predicted by each electrodynamic theory. Alternative experiments that investigate the continuous deformation of fluids subjected to optical radiation pressure may benefit from multi-physics simulation based on statistical mechanics, which are locally non-linear and rigorously conserve momentum. This thesis applies the lattice-Boltzmann method to describe both the electrodynamic and kinetic sub-systems in a multi-physics simulation of a quasi-stationary multiphase fluid subjected to optical radiation pressure. The simulations are used to model multiple electrodynamic formalisms, which yield distinct predictions of microscale fluid deformation. Multiphysics simulations based in statistical mechanics may assist in realizing a single theory for electrodynamics.
Simulations of optical radiationpressure on deformable mediabyCael WarnerB.ASc., The University of British Columbia, 2019A THESIS SUBMITTED IN PARTIAL FULFILLMENT OFTHE REQUIREMENTS FOR THE DEGREE OFMASTER OF APPLIED SCIENCEinTHE COLLEGE OF GRADUATE STUDIES(Electrical Engineering)THE UNIVERSITY OF BRITISH COLUMBIA(Okanagan)December 2021c© Cael Warner, 2021The following individuals certify that they have read, and recommendto the College of Graduate Studies for acceptance, a thesis/dissertation en-titled:Simulations of optical radiation pressure on deformable mediasubmitted by Cael Warner in partial fulfilment of the requirements of thedegree of Master of Applied ScienceKenneth J. Chau, School of EngineeringSupervisorLo¨ıc Markley, School of EngineeringCo-SupervisorJoshua Brinkerhoff, School of EngineeringSupervisory Committee MemberSunny Li, School of EngineeringUniversity ExaminerAdditional Committee Members include:Jonathan Holzman, School of EngineeringSupervisory Committee MemberiiAbstractMultiple electrodynamic theories have existed for over a century, whichhas led to controversy over a field-based description of light-matter inter-actions. Experiments that aim to support or disprove an electrodynamictheory benefit from accurate numerical simulation of each theory’s predic-tions of light-matter interactions. However, most experiments and their nu-merical simulations are limited to cases of rigid-body translation or elasticdeformation, which are equivalently predicted by each electrodynamic the-ory. Alternative experiments that investigate the continuous deformation offluids subjected to optical radiation pressure may benefit from multi-physicssimulation based on statistical mechanics, which are locally non-linear andrigorously conserve momentum. This thesis applies the lattice-Boltzmannmethod to describe both the electrodynamic and kinetic sub-systems in amulti-physics simulation of a quasi-stationary multiphase fluid subjected tooptical radiation pressure. The simulations are used to model multiple elec-trodynamic formalisms, which yield distinct predictions of microscale fluiddeformation. Multiphysics simulations based in statistical mechanics mayassist in realizing a single theory for electrodynamics.iiiLay SummaryThis thesis begins with a review of theoretical and experimental investi-gations into light-matter interactions, and the development of analytical andnumerical approaches for their prediction. Classical light can be described asa population of interacting waves, while classical matter can be described asa population of interacting particles. This thesis applies a discrete numericalmethod for describing the most likely interaction between these populationsof light waves and mass particles. In particular, lattice-Boltzmann numeri-cal methods which invoke Boltzmann’s equation from statistical mechanics.Lattice-Boltzmann numerical methods are used to fulfill the key goals andcontributions of this thesis, which include:1. Comparing the accuracy and computational time of the Yee-latticefinite-difference method and the electrodynamic lattice-Boltzmann methodin predicting a force caused by light;2. Implementing a multi-scale simulation of interactions between lightand a fluid in two and three dimensions; and3. Simulating several varying theoretical predictions for how light canmove or reshape fluids in two and three dimensions.While this thesis attempts to develop a general approach for predicting theinteraction of light and matter, it only provides an initial framework. Fur-ther adaptation is required to describe a specific fluid used in a real-world ex-periment. Nonetheless, the framework is simple to implement and adapt fordifferent experiments given adequate computational resources. The numer-ical method or approach used in this thesis illustrates cases where multipletheories predict either similar or distinct light-matter interactions. Predict-ing light-matter interactions on deformable media may help guide experi-ments that aim to identify a single theory for light-matter interactions. Pre-dicting light-matter interactions is also important for a number of emergingtechnologies which use optical trapping and optical acceleration for micro-assembly or micro-assays.ivPrefaceChapters 3, 4, and 5 include paraphrased information from co-authoredmanuscripts with intended submission to peer-reviewed journals. Through-out this thesis, theories, analytical and numerical methods developed bynumerous authors are utilized. The original authors are appropriately citedwhere their methods are used or described. The combination of these meth-ods and their use for a particular dynamic system are original to this work.vTable of ContentsAbstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iiiLay Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ivPreface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vTable of Contents . . . . . . . . . . . . . . . . . . . . . . . . . . . viList of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xList of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xiList of Acronyms . . . . . . . . . . . . . . . . . . . . . . . . . . . xixAcknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . xxiDedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .xxiiChapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . 11.1 Brief history of optical radiation pressure . . . . . . . . . . . 11.2 Simulations of optical radiation pressure . . . . . . . . . . . . 51.3 The goals and contributions of this thesis . . . . . . . . . . . 61.4 The organization of this thesis . . . . . . . . . . . . . . . . . 6Chapter 2: Theory . . . . . . . . . . . . . . . . . . . . . . . . . . 82.1 Multi-physics energy and momentum conservation . . . . . . 82.2 Electrodynamic sub-system theory . . . . . . . . . . . . . . . 122.3 Kinetic sub-system theory . . . . . . . . . . . . . . . . . . . . 162.3.1 Mass and momentum conservation . . . . . . . . . . . 182.4 Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . 192.4.1 The central difference equation . . . . . . . . . . . . . 202.4.2 Lattice stencils . . . . . . . . . . . . . . . . . . . . . . 20viTABLE OF CONTENTS2.4.3 Discrete spatial operators for the lattice-Boltzmannequation . . . . . . . . . . . . . . . . . . . . . . . . . . 232.4.4 The lattice-Boltzmann equation . . . . . . . . . . . . . 242.4.5 Collision and propagation operators in the lattice-Boltzmannequation . . . . . . . . . . . . . . . . . . . . . . . . . . 25Chapter 3: Comparing finite-difference methods and lattice-Boltzmann methods solutions of kinetic force den-sity . . . . . . . . . . . . . . . . . . . . . . . . . . . . 273.1 Initial condition . . . . . . . . . . . . . . . . . . . . . . . . . . 273.1.1 Hermite Gaussian wave-functions . . . . . . . . . . . . 283.2 The Yee-lattice finite-difference time-domain method . . . . . 313.2.1 Electric and magnetic fields . . . . . . . . . . . . . . . 313.2.2 Energy density and momentum density . . . . . . . . 333.2.3 Force density . . . . . . . . . . . . . . . . . . . . . . . 343.2.4 Boundary conditions . . . . . . . . . . . . . . . . . . . 373.2.5 Limitations . . . . . . . . . . . . . . . . . . . . . . . . 373.3 The electrodynamic lattice-Boltzmann method . . . . . . . . 403.3.1 Electric and magnetic fields . . . . . . . . . . . . . . . 403.3.2 Two-dimensional lattice-matrices . . . . . . . . . . . . 423.3.3 Three-dimensional lattice-matrices . . . . . . . . . . . 443.3.4 Ponderable media . . . . . . . . . . . . . . . . . . . . 443.3.5 Energy density and momentum density . . . . . . . . 453.3.6 Force density . . . . . . . . . . . . . . . . . . . . . . . 463.3.7 Dimensional conversions . . . . . . . . . . . . . . . . . 473.3.8 Limitations . . . . . . . . . . . . . . . . . . . . . . . . 473.4 Comparative simulation configuration . . . . . . . . . . . . . 483.5 Comparative analysis . . . . . . . . . . . . . . . . . . . . . . . 503.6 Computational time and accuracy comparison . . . . . . . . . 503.6.1 Time-averaged force density . . . . . . . . . . . . . . . 53Chapter 4: Multi-physics simulation of optical radiation pres-sure on a multi-phase fluid . . . . . . . . . . . . . . 544.1 Review of multi-physics simulation methods . . . . . . . . . . 544.2 The pseudopotential lattice-Boltzmann method . . . . . . . . 564.2.1 Initial condition . . . . . . . . . . . . . . . . . . . . . 564.2.2 Mass and momentum . . . . . . . . . . . . . . . . . . 584.2.3 The pseudopotential force density . . . . . . . . . . . 584.2.4 The kinetic force density . . . . . . . . . . . . . . . . . 604.2.5 Collision operation and operators . . . . . . . . . . . . 61viiTABLE OF CONTENTS4.2.6 Equilibrium probability mass function . . . . . . . . . 654.2.7 Propagation operation . . . . . . . . . . . . . . . . . . 664.2.8 Boundary conditions . . . . . . . . . . . . . . . . . . . 674.2.9 Dimensionless numbers . . . . . . . . . . . . . . . . . 704.2.10 Dimensional conversion . . . . . . . . . . . . . . . . . 714.3 Algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72Chapter 5: Two- and three-dimensional optical accelerationand deformation of multi-phase fluids . . . . . . . 735.1 Two-dimensional simulation . . . . . . . . . . . . . . . . . . . 745.1.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . 745.1.2 Flat liquid-vapor interface . . . . . . . . . . . . . . . . 775.1.3 Circular liquid-vapour interface . . . . . . . . . . . . . 795.1.4 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 805.2 Three-dimensional simulation . . . . . . . . . . . . . . . . . . 815.2.1 Methodology . . . . . . . . . . . . . . . . . . . . . . . 815.2.2 The linearly-polarized Hermite-Gaussian beam . . . . 835.2.3 The circularly-polarized Hermite-Gaussian beam . . . 855.2.4 Kinetic momentum transfer and electric field polar-ization . . . . . . . . . . . . . . . . . . . . . . . . . . . 875.2.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 89Chapter 6: Conclusion . . . . . . . . . . . . . . . . . . . . . . . . 906.1 Recommendations . . . . . . . . . . . . . . . . . . . . . . . . 93Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94Appendices . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108Appendix A: Analytical methods . . . . . . . . . . . . . . . . . . . 108A.1 Coordinate systems . . . . . . . . . . . . . . . . . . . . . . . . 108A.2 Matrix operations . . . . . . . . . . . . . . . . . . . . . . . . 109A.3 Polynomials and functions . . . . . . . . . . . . . . . . . . . . 111Appendix B: Momentum flux density tensors expanded . . . . . . . 113B.1 Abraham’s momentum flux density tensor . . . . . . . . . . . 113B.2 Minkowski’s momentum flux density tensor . . . . . . . . . . 113B.3 Einstein & Laub’s momentum flux density tensor . . . . . . . 114B.4 Chu’s momentum flux density tensor . . . . . . . . . . . . . . 114B.5 Ampere’s momentum flux density tensor . . . . . . . . . . . . 114B.6 Helmholtz momentum flux density tensor . . . . . . . . . . . 115B.7 Divergence of a stress-tensor . . . . . . . . . . . . . . . . . . . 115viiiTABLE OF CONTENTSB.7.1 Two-dimensional Minkowski stress-tensor divergence . 116B.7.2 Two-dimensional Chu stress-tensor divergence . . . . . 116Appendix C: Dimensional Conversions . . . . . . . . . . . . . . . . 117Appendix D: Program Implementations . . . . . . . . . . . . . . . 118D.1 Electrodynamic lattice-Boltzmann method force density . . . 118D.2 Shan-Chen force density with applicable equation of state . . 120D.3 Periodic boundary conditions . . . . . . . . . . . . . . . . . . 121D.4 Collision operation . . . . . . . . . . . . . . . . . . . . . . . . 121D.5 Propagation operation . . . . . . . . . . . . . . . . . . . . . . 123ixList of TablesTable 2.1 Minkowski formulation of electrodynamics . . . . . . . 12Table 2.2 Einstein & Laub formulation of electrodynamics . . . . 13Table 2.3 The E-H formulation of electrodynamics . . . . . . . . 13Table 2.4 Electrodynamic Momentum and Power Density. . . . . 14Table 2.5 Electrodynamic Momentum Flux Density Tensors. . . 15Table 2.6 Moments of the particle population, n(x,v, t). . . . . . 16Table 2.7 Properties of lattice-stencils used in this thesis . . . . 22Table 3.1 Conversion of dimensionless electrodynamic sub-systemvariables . . . . . . . . . . . . . . . . . . . . . . . . . . 47Table 4.1 Conversion of dimensionless kinetic sub-system variables 71Table A.1 Hermite polynomials . . . . . . . . . . . . . . . . . . . 112Table C.1 Dimensional conversion of system variables . . . . . . 117xList of FiguresFigure 2.1 Electrodynamics of rigid bodies . . . . . . . . . . . . 11Figure 2.2 Illustration of lattice stencils . . . . . . . . . . . . . . 21Figure 2.3 The collision process . . . . . . . . . . . . . . . . . . . 26Figure 3.1 The two-dimensional Yee-lattice . . . . . . . . . . . . 32Figure 3.2 Visual comparison of YL-FDTDM and ELBM numer-ical dispersion . . . . . . . . . . . . . . . . . . . . . . 39Figure 3.3 Discrete lattice-vectors representing a p-polarized elec-tric field . . . . . . . . . . . . . . . . . . . . . . . . . 42Figure 3.4 Discrete lattice-vectors representing an s-polarized elec-tric field . . . . . . . . . . . . . . . . . . . . . . . . . 43Figure 3.5 Illustration of the spatial domain comparing Yee-latticeand lattice-Boltzmann methods . . . . . . . . . . . . 49Figure 3.6 Computational time of the YL-FDTDM and the ELBM 51Figure 3.7 Numerical accuracy of the YL-FDTDM and the ELBM 52Figure 3.8 Ratios of errors and computational times with respectto spatial units per wavelength . . . . . . . . . . . . . 52Figure 3.9 Chu time-averaged force density in the YL-FDTDMand the ELBM . . . . . . . . . . . . . . . . . . . . . . 53Figure 4.1 Liquid-vapor phase coexistence . . . . . . . . . . . . . 60Figure 4.2 Illustration of periodic boundary conditions . . . . . . 67Figure 4.3 Bounce-back boundary conditions . . . . . . . . . . . 68Figure 4.4 Non-equilibrium bounce-back boundary conditions . . 70Figure 4.5 Multi-physics simulation algorithm . . . . . . . . . . 72Figure 5.1 Two-dimensional simulation configuration used to ex-amine the deformation of a flat liquid-vapour interface 75Figure 5.2 Two-dimensional simulation configuration used to ex-amine the deformation of a circular liquid-vapour in-terface . . . . . . . . . . . . . . . . . . . . . . . . . . 76xiLIST OF FIGURESFigure 5.3 Two-dimensional dielectric liquid-vapor interface de-formation from a p-polarized electric field . . . . . . . 77Figure 5.4 Two-dimensional dielectric liquid-vapor interface de-formation from an s-polarized electric field . . . . . . 78Figure 5.5 Two-dimensional dielectric droplet deformation . . . 79Figure 5.6 Two-dimensional magnetic droplet deformation . . . 80Figure 5.7 Three-dimensional simulation configuration used toexamine the deformation of a spherical liquid-vapourinterface . . . . . . . . . . . . . . . . . . . . . . . . . 82Figure 5.8 Three-dimensional deformation of an incompressiblespherical droplet from a linearly-polarized Gaussianbeam . . . . . . . . . . . . . . . . . . . . . . . . . . . 84Figure 5.9 Three-dimensional deformation of an incompressiblespherical droplet from a circularly-polarized Gaussianbeam . . . . . . . . . . . . . . . . . . . . . . . . . . . 86Figure 5.10 Polarization-dependent momentum transfer accord-ing to the Chu and Amperian/Lorentz formalisms . . 87Figure 5.11 Polarization-independent momentum transfer accord-ing to the Abraham, Minkowski, Einstein & Laub,and Helmholtz formalisms . . . . . . . . . . . . . . . 88Figure A.1 Difference between the position vector and displace-ment vector . . . . . . . . . . . . . . . . . . . . . . . 109xiiList of Symbols(x0, y0, z0) Position of the beam-waist on the beam-axis, page 29[n], [nˆ] Population density as a vector in velocity space and mo-ment space, page 63[Λ], [Λˆ] Collision matrix in velocity space and moment space,page 63[M] Transformation matrix, page 63ℵ Refractive index, page 28α, β, γ, µ, ν Order and dimensional indices, page 18G Momentum density, page 8S Power per cross-sectional area, page 8∆t Dimensional time unit, page 36δt Dimensionless time unit, page 25∆x Dimensional spatial unit, page 33δx Dimensionless spatial unit, page 25 Permittivity of a non-vacuum medium, page 140 Free-space permittivity, page 12r Relative permittivity, page 14η Impedance of the source medium, page 30nˆi Normal unit vector at an interface, page 57ι Interface thickness, page 57λ0 Free-space wavelength, page 28xiiiLIST OF SYMBOLSR2 Determination coefficient, page 79k Wave vector, page 11T Torque, page 11B Magnetic flux density, page 12D Electric flux density, page 12E Electric field intensity, page 12e Discrete electric field lattice vector, page 41f Force density, page 8fSC Shan-Chen force density, page 59H Magnetic field intensity, page 12h Discrete magnetic field lattice vector, page 41i Imaginary unit, page 29JE Electric current density, page 12JM Magnetic current density, page 12Meq Magnetization density equilibrium state, page 44Peq Polarization density equilibrium state, page 44r Displacement vector, page 57ri Locus of the interface with respect to the origin, page 57s Two-dimensional axis vector for a Hermite-Gaussian beam,page 29s⊥ Two-dimensional curvature vector for a Hermite-Gaussianbeam, page 29u Mean velocity of a particle population, page 13v Random velocity of a particle population with referenceto a stationary absolute space, page 16x Vector position, page 16xivLIST OF SYMBOLSR Gas constant, page 16T Mean temperature field, page 16µ Permeability of a non-vacuum medium, page 14µ0 Free-space permeability, page 12µr Relative permeability, page 14ν Kinematic viscosity, page 71 Into page, page 43Ω Collision operator, page 17ω Angular frequency, page 28⊗ Out of page, page 43T momentum flux density tensor, page 8Φ Power density, page 8φnm Hermite-Gaussian wave-function, page 28ρ Mean mass density, page 16ρc Two-phase critical point mass density, page 60ρv Electric charge density, page 12ρinit. Initial mass density, page 57ρliq. Liquid mass density, page 57ρvap. Vapour mass density, page 57σE Electric conductivity, page 12σM Magnetic conductivity, page 12σαβ Stress tensor, page 19σk,xx, σk,xy Diagonal and off-diagonal components of the Cauchy stress-tensor in two dimensions, page 64τ Discrete relaxation time, page 25xvLIST OF SYMBOLSMa Mach number, page 77pk,x Net dimensionless kinetic momentum in the xˆ direction,page 87pk,y Net dimensionless kinetic momentum in the yˆ direction,page 87pk,z Net dimensionless kinetic momentum in the zˆ direction,page 87θ(z) Divergence angle of the Hermite-Gaussian beam, page 29ϑ Angle between the beam-axis and the source-axis, page 29v˜ Random velocity of a particle population relative to itsmean velocity, page 16ζ Bulk viscosity, page 63ᵀ The transpose operation, page 64a, b, R Constants used in the Carnahan-Starling equation of state,page 59c Speed of light in a vacuum, page 28c0 Dimensionless speed of light, page 45cs The lattice speed of sound, page 25e Local energy density per mass density, page 16E0 Peak electric field intensity, page 30f Generic function, page 20G Interaction potential (constant), page 59gm Electrodynamic psuedoparticle probability mass function,page 40g∗m Post-collision pseudo-particle population, page 40geqm Equilibrium electrodynamic psuedoparticle probability massfunction, page 40xviLIST OF SYMBOLSi Spatial index representing a column element in a discretefield, page 23j Spatial index representing a row element in a discretefield, page 23k Wave-number, page 28kB Boltzmann constant, page 17l ∈ {s, k, e} Index for thermodynamic (s), kinetic (k), or electrody-namic (e) sub-system, page 8mp Particle mass, page 17n Particle population density, time-index, or modal index,page 16neq Equilibrium particle population density, page 16nλ Number of spatial units per wavelength, page 50Nm Number of collision outcomes for the particle populationprobability mass function, page 24nm Particle population as a probability mass function, page 24n∗m Post-collision particle population, page 25neqm Equilibrium particle population as a probability massfunction, page 24NT Number of time-periods, page 36Nt Number of time steps, page 36Nx Number of spatial lattice units in the xˆ-direction, page 39Ny Number of spatial lattice units in the haty-direction,page 39p Pressure derived from an equation of state, page 59q(z) Complex beam parameter, page 29R(z) Beam radius of curvature, page 29xviiLIST OF SYMBOLSs Polar axis coordinate, page 29Sm Forcing term used in the collision operation, page 61sm Relaxation parameters in a multi-relaxation approach,page 63T Time period, page 28t Time, page 16V Generic variable representing a scalar field, page 12W Energy density, page 8w(z) Beam-waist of Hermite-Gaussian beam, page 29we,m Lattice weights for an electrodynamic sub-system, page 22wk,m Lattice weights for a kinetic sub-system, page 22x Dimension that is horizontal on the page in this thesis,page 23y Dimension that is vertical on the page in this thesis,page 23z Out-of-plane dimension (2D), vertical dimension (3D),page 29zR Rayleigh range, page 29d Number of dimensions, page 16xviiiList of AcronymsAB Abraham. 14, 15, 77–80, 83–86, 88, 89, 92AMP Amperian. 14, 15, 77–80, 83–87, 89BE Boltzmann’s equation. 17BGK Bhatnagar-Gross-Krook. 25, 61, 65, 66, 71BPML Berenger perfectly matched layer. 37, 48CET Chapman-Enskog theory. 19CHU Chu. 14, 15, 77–80, 83–87, 89CS-EOS Carnahan Starling equation of state. 60D1Q3 one-dimensional three lattice vector. 19, 22,23, 46D2Q5 two-dimensional five lattice vector. 19–23, 42,43, 46D2Q9 two-dimensional nine lattice vector. 19, 21–23, 25, 26, 62–64, 69, 71, 75, 76D3Q27 three-dimensional twenty-seven lattice vector.19, 21–23, 25, 71D3Q7 three-dimensional seven lattice vector. 19–23,40, 46EL Einstein & Laub. 14, 15, 77–80, 83–86, 88, 89ELBM electrodynamic lattice-Boltzmann method.xi, 7, 19, 27, 38–40, 44–54, 56, 72, 73, 81,118EOS equation of state. 16, 19EPMF equilibrium probability mass function. 56, 61,65, 66FDTD finite-difference time-domain. xxi, 5, 54xixList of AcronymsFEM finite-element method. 5HLM Helmholtz. 14, 15, 77–80, 83–86, 88, 89KBC Karlin-Boltzmann Collision. 65, 66LBE lattice-Boltzmann equation. 19, 20, 24, 66LBM lattice-Boltzmann method. 20, 55, 59, 65–67,70MDW mass-density wave. 4MN Minkowski. xii, 13–15, 77–80, 83–86, 88, 89,92, 116MRT multi-relaxation time. 62, 64–66, 76NEBB non-equilibrium bounce-back. 75NSE Navier-Stokes Equation. 55, 66PBC periodic boundary conditions. 82PLBM pseudopotential lattice-Boltzmann method.7, 19, 55, 56, 59, 61, 65, 66, 70, 73–81, 118PMF probability mass function. 56, 58, 59, 61, 63,66, 67SEM stress-energy-momentum. 4, 8, 15SRT single-relaxation time. 64–66, 75, 76, 81TE transverse electric, s-polarized. 32, 43TEM transverse Hermite-Gaussian. 28–30TM transverse magnetic, p-polarized. 32, 35, 42YL-FDTDM Yee-lattice finite-difference time-domainmethod. xi, 7, 19, 27, 31, 34, 36–40, 44,47–53, 72, 73, 81xxAcknowledgementsI would like to extend my deep appreciation to the Natural Sciencesand Engineering Research Council (NSERC) of Canada and the Universityof British Columbia (UBC) College of Graduate Studies for their generoussupport in the form of the NSERC Canadian Graduate Scholarship - Masters(CGSM), the NSERC Discovery Grant (DG), and the UBC Graduate DeansEntrance Scholarship (GDES). This work was also made possible with theassistance of Dr. Lo¨ıc Markley and Dr. Kenneth J. Chau, who were availablefor bi-weekly discussions on the topic of optical radiation pressure through-out the COVID-19 pandemic. Lo¨ıc provided helpful suggestions for theelectrodynamic sub-system, including improvements to the initial simula-tion methodology that increased its performance. He also brought a healthyskepticism to the methodology that helped it mature. Kenneth, with theassistance of his former student Max Bethune-Waddell, also provided theinitial framework for computational electrodynamics that has enabled thisthesis. In particular, the finite-difference time-domain (FDTD) methods ofsolving kinetic force density, and extensive testing of the force density ofvarious electrodynamic formalisms on rigid and deformable media. I wouldfurther like to thank Brock Legentil for his time as an undergraduate re-search assistant, whose quick learning, impressive attention to detail andcomputational resources assisted our simulations of two-dimensional FDTDforce-densities. Andreas Hauser and Dr. Chun-Sheng Wang, thank-you foryour guidance and suggested practices with the electrodynamic and multi-phase lattice-Boltzmann methods. Special thanks to all anonymous review-ers and editors who helped clarify this work. Finally, I would like to thankVincent Loi for showing me how to script acronymns and symbols in thethesis template.xxiDedicationI dedicate this work to all researchers who endeavoured to investigatelight-matter interactions; your theories have stood the test of time. In par-ticular, James C. Maxwell, who formulated much of the fundamental theoryfor both the electrodynamic and kinetic sub-systems explored in this thesis.Among Maxwell were Ludwig Boltzmann, William Thomson (Lord Kelvin),Hendrik Lorentz, John Henry Poynting, Albert Einstein, Jakob Laub, MaxAbraham, Hermann Minkowski, Lan Jen Chu, Robert Fano, and Hermannvon Helmholtz among many other scientists whose theories fundamentallycontribute to this thesis. Further, I dedicate this work to those who haveprovided support to researchers who investigate optical radiation pressure;thank-you for appreciating the significance, challenges and value of under-standing these ubiquitous forces acting at a distance.More personally, I dedicate this thesis to Alice Jaquish for her kind wordsof encouragement as I pursued a university education; you may be gone, butyou will be remembered. I also dedicate this thesis to Jan Cyganik; he wasa good man who cared much for the well-being of others. I further dedicatethis thesis to Josh, Brandon, Aspen and Joe; thank-you for the good times.xxiiChapter 1IntroductionThis chapter discusses the brief history of optical radiation pressure, fromits discovery to its present theoretical debate. This brief history is followedby the computational methods used to predict optical radiation pressure,and how such methods are used to understand the phenomenon. Last, thischapter discusses the organization of subsequent chapters in this thesis witha short description of their content.1.1 Brief history of optical radiation pressureThe first recorded observation of optical radiation pressure wasmade by James Kepler in 1679, when he observed a comet’s tail directedaway from the sun [1, 2]. It was almost two centuries later, in 1845, whenelectrostatic and magnetostatic forces on ponderable media were predictedand measured by Lord Kelvin [3]. In 1865, a dynamic relationship betweenelectric and magnetic fields was discovered by Michael Faraday and devel-oped into a classical theory of electromagnetism by James C. Maxwell [4].James C. Maxwell’s classical theory of electromagnetism describes the cou-pling of electric and magnetic fields as electromagnetic waves that manifestas light. Maxwell’s equations imply that light can exert force on matterand therefore, it lets us conclude that light carries momentum. In 1881,Helmholtz defined electric and magnetic forces that satisfy Lord Kelvin’smeasurements, and in 1895, Lorentz derived electric and magnetic forcesfrom Maxwell’s equations [5, 6]. Although Helmholtz and Lorentz had in-troduced unique expressions for an electromagnetic force density, the formerbased on a polarization density and the latter based on free-charge, theseforce density expressions describe the same net force acting on a ponderablemedium [7, 8]. In 1901, Lebedew confirmed the existence of optical radi-ation pressure by measuring the displacement of a solid dielectric [9]. In1905, Poynting proposed that optical radiation pressure acts normal to aboundary [10].11.1. Brief history of optical radiation pressureMultiple electrodynamic formalisms were developed shortly afterthe theory of special relativity. The theory of special relativity, which wasdeveloped by Albert Einstein in 1905, suggests that mass, length, and timeare velocity-dependent variables that change with respect to a universallyconstant vacuum speed of light [11]. Lorentz discovered that Maxwell’sequations also satisfy special relativity, given the appropriate Lorentz trans-formations for electric and magnetic fields in an arbitrary frame of reference[12]. Special relativity established the need for a frame-invariant, or Lorentzinvariant, electrodynamic formalism capable of conserving energy and mo-mentum in all frames of reference. In 1908 and 1909, Hermann Minkowski,Max Abraham, and Albert Einstein & Jakob Laub derived unique electro-dynamic formalisms from the relativistic Maxwell’s equations, energy andmomentum continuity [13–15]. Each formalism satisfies energy and momen-tum continuity, and therefore predicts the same net optical force on a rigidbody [7, 8]. Minkowski and Abraham’s formalisms also describe the Lorentzforce density, despite their distinct forms of momentum [13, 14]. Einstein& Laub presented a unique force density in their formalism, such that thetorque density includes light-matter interactions exclusive of the force den-sity [15, 16].Early experiments with optical radiation pressure made definitiveconclusions in favour of a single electrodynamic formalism that were latercriticized and found inconclusive [7, 8, 17]. Theoretical controversy inten-sified given that each electrodynamic formalism has been supported in ex-periment [7, 17, 18]. Many experiments analyzed the deformation of elasticsolids such as cantilevers, torsion beams, spectral mirrors, and optical fibers[17]. While the conclusions of these studies argue in favour of a single electro-dynamic theory, the net forces and kinematics measured in each study havebeen identically predicted by each electrodynamic formalism in numericalanalyses [7, 8].Alternative electrodynamic theories were developed in the 1960’sin an attempt to resolve the ensuing controversy that had lasted for overhalf a century [5, 19, 20]. Lan Jen Chu developed the E-H formulation ofMaxwell’s equations in the 1960’s, which is a microscopic formulation thatdepends only on external electric and magnetic field intensity [20]. Usingthe E-H formulation, with principles of energy and momentum conservation,Chu developed an electrodynamic formalism that describes the Lorentz forcedensity and Abraham momentum density from electric and magnetic field21.1. Brief history of optical radiation pressureintensities alone. Since Chu’s formalism does not explicitly include the elec-tric and magnetic dipoles of ponderable media, it is philosophically distinctfrom the Abraham, Minkowski, and Einstein & Laub formalisms. Alter-natively, the Lorentz/Amperian formalism described by Robert Fano onlydiffers from Chu’s formalism in describing magnetic dipoles in magnetic me-dia [19].Hidden momentum became a topic of the controversy following itsdiscovery by Shockley in 1967 [21]. Mansuripur demonstrated that the Chuand Lorentz/Amperian formalisms for a dielectric neglect the influence of“hidden momentum,” in favour of the Lorentz force density [16, 21]. Conse-quently, there may exist additional terms for light-matter interactions in anelectromagnetic torque density that are exclusive of an electromagnetic forcedensity [16]. Although Einstein may have considered these additional termserroneous in the formalism he developed with Laub, Mansuripur demon-strates that these additional terms may account for the “hidden momentum”discovered by Shockley in 1967 [16, 21]. Plotting the Lorentz force densitythroughout a dielectric medium, Mansuripur also predicted that it may causedeformation opposite to earlier experiments on fluids performed by Ashkin& Dziedzic [22]. Such a conclusion was also drawn by Bethune-Waddell &Chau who demonstrated that the Lorentz/Amperian and Chu formalisms, inaddition to predicting electric-field polarization-dependent optical radiationpressure, can predict opposite interface deformation and optical trappingcompared with experiment [7, 8]. Given that the re-interpreted Einstein& Laub formalism may accurately predict the deformation of generic fluidssubjected to optical radiation pressure, its predictions for the electrodynam-ics of elastic media have also been explored [8, 16, 23].Modern experiments with optical radiation pressure investigateits deformation of microscopic fluids or nanometric solids. Optical radiationpressure on ponderable fluids was first studied by Ashkin & Dziedzic in 1973,who subjected an air-water interface to a focused Gaussian laser beam [22].Highly deformable toluene-water interfaces have also been investigated, andsuggested as a means of discriminating between the predictions of electrody-namic formalisms [8, 24–27]. Ashkin & Dziedzic’s experiment was repeatedin 2014 by Astrath et al., with improved measurement sensitivity and preci-sion, and demonstrated agreement between the measured deformation andthe simulated prediction of the Helmholtz formalism [5, 28]. However, As-trath et al.’s study did not compare measurements of the fluid-fluid inter-31.1. Brief history of optical radiation pressureface deformation with the predicted deformation of other electrodynamicformalisms [23]. Zhang et al. demonstrated that fluid depth is also animportant consideration in experiments deducing the Minkowski or Abra-ham momenta, since the trivial solution of the Laplace pressure becomesdominant in the Navier-Stokes equation with decreasing fluid depth [29].However, Zhang et al.’s study approximates pressure at an interface accord-ing to Abraham and Minkowski momenta alone, while there exist alternativeformalisms that were not compared and that may provide an alternative ex-planation for the downward interface deflection [7]. Pozar et al. comparedpredictions of acoustic vibrations from multiple electrodynamic formalismswith their piezoelectric measurement from pulsed optical radiation pressureon an elastic dielectric mirror [23]. Comparisons between the measured andsimulated acoustic waves demonstrate that Chu and Amperian formalismspredict significantly less vibration than had been measured [23]. However,Pozar et al. noted that increased measurement precision may be requiredto discriminate between the remaining formalisms [23].The theoretical foundation of electrodynamics remains contested,especially in studies where the frame invariance of multiple electrodynamicformalisms are re-examined. While each electrodynamic formalism was orig-inally derived to satisfy Lorentzian and Lagrangian invariance, to ensuretheir independence from a particular frame of reference, their invariance withthe frame of reference remains a subject of controversy [30, 31]. Partanenand Tulkki performed experiments of radiation pressure on thin films anddeveloped a mass-polariton theory which remedies the Lagrangian variancein Minkowski’s stress-energy-momentum (SEM) tensor [30]. In the mass-polariton theory, mass-polaritons are quasi-particles which travel at thespeed of sound within a medium, and represent interactive coupling betweenphotons and atoms. Mass-polaritons produce an atomic mass-density wave(MDW) with acoustics described by Minkowski’s momentum [30]. Includingthe MDW and the rest energy of a material system in the energy continu-ity equation modifies the Minkowski SEM tensor to satisfy Lorentzian andLagrangian invariance [30]. The mass-polariton theory also explains the dif-ference in center-of-mass translation predicted by Minkowski and Abrahammomenta, and how Minkowski’s momentum density satisfies Hamiltonianmechanics [32, 33]. In response, Kemp et al. emphasized the non-causalityof the Minkowski SEM tensor [31]. Sheppard & Kemp also argue that theAbraham, Einstein & Laub, and Lorentz/Amperian SEM tensors are La-grangian variant, such that Chu’s SEM tensor is the only formalism which41.2. Simulations of optical radiation pressureis Lagrangian invariant [31]. Therefore, the theoretical foundations of opticalradiation pressure remain a topic of debate since there exist multiple theo-retical interpretations of each electrodynamic formalism’s frame-invariance,energy and momentum conservation [34]. Resolving the debate will requireeither increased precision and accuracy measurement or more experimentselucidating distinct predictions from each of the formalisms [23]. Such ex-periments may benefit from the simulated deformation of quasi-stationarymedia in a laboratory frame of reference [7, 8, 23, 28].1.2 Simulations of optical radiation pressureSimulations of optical radiation pressure using finite-difference time-domain (FDTD) methods have been successful in representing multiple elec-trodynamic experiments on rigid media [7, 8, 33, 35]. FDTD has also beencombined with front-tracking or the finite-element method (FEM) to sim-ulate the deformation of liquids and viscoelastic media in optical radiationpressure experiments [7, 8, 23, 28]. Alternative numerical methods includethe boundary-integral-element method with ray tracing to approximate op-tical radiation pressure [25]. Recently, there has been increased interest inlattice-Boltzmann methods for electrodynamics and multi-phase fluids [36–40]. Not only are electrodynamic and kinetic lattice-Boltzmann methodssimpler to implement than FDTD or solvers of Maxwell’s equations or theNavier-Stokes equation, but they can more easily represent natural inter-faces and mesoscopic phenomena [41]. Mesoscopic phenomena occur at amesoscopic scale of condensed matter physics. The mesoscopic scale liesbetween the microscopic scale (which considers individual interacting atomsand molecules) and the macroscopic scale (which considers continuous ma-terials). In addition, the kinetic lattice-Boltzmann method can simulatenon-linear hydrodynamics locally, by-passing the need for iterative solutionsof hyperbolic differential equations that are assumed linear in finite-elementmethods [41–43]. While inefficient for modeling steady flow, kinetic lattice-Boltzmann methods are more efficient for modeling weakly compressibleflows [44]. The lattice-Boltzmann method has also been adapted to representacoustics and high-Reynolds number, compressible, thermal and relativisticfluid flows [45–50].51.3. The goals and contributions of this thesis1.3 The goals and contributions of this thesisThe goals and contributions of this thesis are to:1. Develop a simple, accurate, and fast electrodynamic simulator of time-averaged optical force density in two and three dimensions.2. Develop a simple, accurate, and re-configurable multi-physics simula-tor of optical radiation pressure on continuously deformable media intwo and three dimensions.3. Simulate optical radiation pressure on a spherical liquid droplet sus-pended in its own vapor to illustrate distinct deformations of a meso-scopic fluid predicted by multiple electrodynamic formalisms.Numerically simulating the mesoscopic optical deformation of ponderablemedia may help guide future experiments aiming to contrast the predic-tions of multiple electrodynamic formalisms [22, 28]. Such simulations maybenefit predictions of micro-assembly with optical tweezers [51], optical ma-nipulation of cells and vesicles for their analysis in non-invasive surgery[52], or other microfluidic assays [36, 38, 53]. With further description ofphotothermal absorption, additional applications include the prediction offluid deformation from plasmonic heating [54], nitrogen bubble nucleationand cavitation in liquid cryogenics [55, 56], or laser focusing [57]. Givenadaptation for arbitrary frames of reference [50, 58, 59], such simulationmethodology may benefit secondary applications in orbital debris removal[60], and interstellar space travel [2].1.4 The organization of this thesisThis thesis uses formal classification to organize its content into chap-ters. The following chapters are organized in order of the thesis goals andcontributions:− Chapter 2, Theory: Describes classical electrodynamic theory, mo-mentum and energy conservation. General conclusions for the elec-trodynamics of rigid bodies, the advantage of a statistical mechanicsapproach for mesoscopic analysis, and general/fundamental numericalmethods and approaches are presented.− Chapter 3, Comparing finite-difference methods and lattice-Boltzmann methods solutions of kinetic force density: The61.4. The organization of this thesisYee-lattice finite-difference time-domain method (YL-FDTDM) and asingle-lattice electrodynamic lattice-Boltzmann method (ELBM) arecompared for their numerical modeling of kinetic force densities innon-dispersive media.− Chapter 4, Multi-physics simulation of optical radiation pres-sure on a multi-phase fluid: The pseudopotential lattice-Boltzmannmethod (PLBM), with a Carnahan-Starling equation of state, is in-tegrated with the ELBM for two-dimensional and three-dimensionalsimulations of multiphase fluid optical acceleration and trapping.− Chapter 5, Two- and three-dimensional optical accelerationand deformation of multi-phase fluids: Simulations of optical ra-diation pressure on mesoscopic liquid droplets suspended in vapor areexplored in both two and three dimensions for multiple electrodynamicformalisms.− Chapter 6, Conclusion: This thesis concludes with a description ofthe general observations from purely kinetic approaches which conservemomentum, and general properties of a quasi-stationary experiment orsimulation thereof, which may benefit present understanding of opticalradiation pressure. Limitations to the study and recommendations forfuture work are emphasized.− Appendix: Contains the relevant mathematical operations and deriva-tions in this thesis and some MATLAB functions for the numericalmethods employed.7Chapter 2TheoryThis chapter describes the interaction between electrodynamic and ki-netic sub-systems that conserve energy and momentum in the rest-frame.This chapter omits descriptions of Lagrangian-invariance or Lorentz-covarianceof SEM tensors, which are outside the scope of this thesis [30, 31]. Maxwell’sequations and six popular electrodynamic formalisms are described in thequasi-stationary or “near-rest” frame in a discussion of the electrodynamicsub-system. Boltzmann’s equation and the discrete lattice-Boltzmann equa-tion are detailed for the kinetic sub-system. Finally, this chapter concludeswith a general description of discrete numerical methods employed in thisthesis. Refer to “analytical methods,” Appendix A, for a discussion of therelevant mathematical operators, polynomials and functions.2.1 Multi-physics energy and momentumconservationIn the rest-frame, a physical system can be separated into three sub-systems indexed by l ∈ {s, k, e}: thermodynamic (l = s), kinetic (l =k), and electrodynamic (l = e) sub-systems. These sub-systems interactaccording to the continuous conservation of their total energy density, W ,and momentum density, G,∇ ·S + ∂∂ tW =∑lΦl, (2.1)∇ ·T + ∂∂ tG =∑lfl, (2.2)where S is the power per cross-sectional area and T is momentum fluxdensity tensor which defines the state of stress at a point for a stationarymedium [8, 20]. In an isolated system, the net power density,∑l Φl, is zeroand the net force density,∑l fl, is a null vector. Energy and momentum are82.1. Multi-physics energy and momentum conservationconserved in their transfer between sub-systems [20],∑lΦl = Φe + Φk + Φs = 0,∑lfl = fe + fk + fs = 0. (2.3)In this thesis, the thermodynamic sub-system is neglected such that Φs = 0and fs = 0 to represent ideal light-matter interactions in an isolated isother-mal system. Modern experiments investigating optical radiation pressure ondeformable liquids are often designed to achieve such conditions, for exam-ple, the case of green laser light transmitted through pure liquid water [28].Consequently, Newton’s third law states that the force density acting on theelectrodynamic sub-system must be equal and opposite to the force densityacting on the kinetic sub-system, such that fk = −fe [7, 8, 20]. The result isa single momentum continuity equation describing the kinetic force densityacting on an isothermal kinetic sub-system in terms of an electromagneticmomentum flux density tensor, Te, and an electromagnetic momentum den-sity, Ge,∇ ·Te + ∂∂ tGe = −fk. (2.4)Multiple theoretical formulations of Te and Ge have been proposed to satisfythis relation, particularly in the rest-frame. To limit comparative studies inthe rest-frame, only a subset of six mathematically debatable frame-invariantelectrodynamic formalisms are compared.In typical optical radiation pressure experiments on macroscopic media,the velocity of the laboratory frame and/or the medium under considerationis negligible in comparison to the speed of light. The medium can then beconsidered quasi-stationary and Lorentz transformation unnecessary [8, 23].The simulations in this thesis assume the medium is quasi-stationary andcontinuously deformable. This is a necessary development for distinguishingbetween electrodynamic formalisms as previous work has shown that theyyield identical rigid-body motion [7, 8]. Previous studies reveal some generalconclusions for the dynamics of magnetodielectric bodies subjected to quasi-stationary optical radiation pressure:1. In one dimension, only linear pressure acting along the direction ofpropagation is accessible. Neither shear stress nor torque can be ana-lyzed in one dimension.2. In one dimension, each electrodynamic formalism predicts distinctforce density but identical time-averaged time-harmonic pressure. Eachformalism therefore predicts identical rigid-body kinematics in mostexperiments [8].92.1. Multi-physics energy and momentum conservation3. In two dimensions, the net line force density acts in the direction ofpropagation [10]. Each formalism predicts unique transverse shearstress resulting from electric and magnetic field attenuation. Thetransverse shear stress can be distinct for each electric field polar-ization.4. Consistent net force on a body consisting of heterogeneous or homo-geneous media can be evaluated when; (i) a source of radiation isexternal to the body, and (ii) the body is fully bounded by free space.Omitting a free-space/medium boundary, and/or initializing a sourceinternal to the body, will result in an inconsistent evaluation of netforce from a distributed force density.5. Given a source of radiation with a symmetric transverse profile, forcesacting orthogonal to the direction of propagation are balanced in eachdirection. Electromagnetic torque is generally non-zero and dependson the electric field polarization and isotropy of the media. Birefrin-gent media will be subjected to electromagnetic torque [61].6. If the body is distributed asymmetrically with respect to an externalsource of radiation, it can be subjected to a net force and/or torquefrom scattering transverse to the direction of incidence.Figure 2.1 illustrates the general conclusions for the dynamics of magnetodi-electric bodies subjected to quasi-stationary optical radiation pressure thatwere described in the previous list.102.1. Multi-physics energy and momentum conservation1 23 45 6zzk zkzk zk=z1-Dimensional1-Dimensional2-Dimensional 2-Dimensional3-Dimensional 3-Dimensionalyk=Τ 0zykΤ Τzkykzkykzzz zxk xk, ,, ,, ,0 0,  0 0,  0 0,  0 0, 0 0,  0 0,  0 0,  0 0, 0 0,  0 0, Figure 2.1: [Color online] Electrodynamics with different source and objectconfigurations in one, two, and three dimensions. In the top left-hand cor-ner of each box is a number representing each general conclusion. Dashedblack lines are major axes, while solid black arrows represent propagationdirections and wave-vectors, {kx,ky,kz}. Solid black lines represent dielec-tric or system boundaries. (1) Neither transverse force density nor shearstress (solid blue arrow) can be evaluated in one dimension. Electrody-namic torque, T, cannot be evaluated in one dimension. (2) Force densityand optical radiation pressure (solid green arrows) can be evaluated in onedimension. The same time-averaged time-harmonic radiation pressure ispredicted by each electrodynamic formalism. (3) Transverse force densityor shear stress can be evaluated in two dimensions. (4) Net line-force densitycan only be evaluated using a source of electromagnetic radiation positionedin free-space (along the green dashed line for example). Net line force densitycannot be evaluated using a source of electromagnetic radiation positionedwithin the medium (along the red dashed line for example). (5) Opticalacceleration and direction of net force (green arrow) in a ponderable spheresubject to optical radiation with a symmetric transverse profile. (6) Netforce and/or torque on a body asymmetrically distributed with respect tothe transverse symmetric profile in (5).112.2. Electrodynamic sub-system theory2.2 Electrodynamic sub-system theoryThe electrodynamic sub-system in a rest frame is typically represented bythe Minkowski formulation of Maxwell’s equations, where the dot notation,i.e. “V˙ ”, is a shorthand for the partial derivative of a field V with respectto time (∂ V/ ∂ t) [4].Table 2.1: Minkowski formulation of electrodynamicsName of Law Vector calculus formGauss’s Law ∇ ·D = ρv∇ ·B = 0Faraday’s Law ∇×E = −B˙− JMAmpere’s Law ∇×H = D˙ + JEIn the Minkowski formulation of Maxwell’s equations, electric and mag-netic field intensities are E and H, while electric and magnetic flux densitiesare D and B, respectively. Dimensional constants 0 and µ0 represent thefree-space permittivity and permeability, respectively. Gauss’s law predictsthe electric charge density, ρv, while free electric and magnetic current den-sitiesJE = σEE and JM = σMH, (2.5)are represented by electric conductivity σE , and magnetic conductivity σM .Although a magnetic current density has yet to be discovered, it servesas a useful analytical tool that satisfies Maxwell’s equations, energy andmomentum conservation [62].Alternative formulations for Maxwell’s equations, which differ philosoph-ically from the Minkowski formulation, lead to alternative predictions of mi-croscopic light-matter interactions. Einstein & Laub’s formulation describesa dependence between E and an electric polarization density P, and H anda magnetic polarization density M, in a magnetodielectric medium. Accord-ing to Einstein & Laub’s formulation, electric and/or magnetic fields alwaysexist in magnetodielectric media, even in the absence of external electricand magnetic fields [15, 16, 63]. In the latter case, the electric and/or mag-netic forces are at equilibrium from the orientation of electric and/or mag-netic dipoles. Einstein & Laub’s formulation (Tab. 2.2) is mathematicallyidentical to Minkowski’s formulation (Tab. 2.1) given constitutive relationsD = 0E + P and B = µ0(H + M). However, Einstein & Laub’s formula-tion establishes a local dependence between field intensity and polarization122.2. Electrodynamic sub-system theorydensity, such that P and M may describe the internal electric and magneticfield intensity distributed throughout a magnetodielectric medium.Table 2.2: Einstein & Laub formulation of electrodynamicsName of Law Vector calculus formGauss’s Law ∇ ·0E = −∇ ·P + ρv∇ ·µ0H = −∇ ·µ0MFaraday’s Law ∇×E = −µ0(H˙ + M˙)− JMAmpere’s Law ∇×H = (0E˙ + P˙) + JEAnother formulation of Maxwell’s equations that considers the generalcase of polarization and magnetization density in fluids moving at relativisticvelocities is called the E-H formulation. Although the E-H formulation ofMaxwell’s equations is consistent with the Minkowski formulation [20], itdirectly includes the velocity of the polarization and magnetization densitywith respect to only external electric and magnetic fields.Table 2.3: The E-H formulation of electrodynamicsName of Law Vector Calculus FormGauss’s Law ∇ ·0E = −∇ ·P + ρv∇ ·µ0H = −∇ ·µ0MFaraday’s Law ∇×E = −µ0(H˙ + M˙)−∇×(µ0M× u)− JMAmpere’s Law ∇×H = (0E˙ + P˙) +∇×(P× u) + JEIn the E-H formulation, E and H are generally different in polarizableor magnetizable media, as compared with their definition in the Minkowskiformulation, which is an important consideration for their Lorentz covariance[20]. The relationship between the field intensities and field densities in thetwo formulations isEMin. = E + µ0M× u, DMin. = 0E + P,HMin. = H−P× u, BMin. = µ0(H + M),(2.6)where u is the mean velocity of particles in the medium, and subscript “Min.”designates the Minkowski formulation [20]. According to Eq. 2.6, local rel-ative velocity of P and M are responsible for variation of EMin. and HMin.from E and H in the free-space Maxwell flux-density tensor. Should themean velocity of particles in a medium have a magnitude that is negligi-ble with respect to the speed of light (|u| c), E and H are the same132.2. Electrodynamic sub-system theoryin both the E-H and Minkowski formulations of Maxwell’s equations [20].This condition is satisfied by the quasi-stationary deformation of media fromtime-averaged time-harmonic electromagnetic waves [8, 23]. Each formula-tion of Maxwell’s equations can be expanded via dyadic products in therest-frame to recover an electrodynamic momentum flux density tensor andmomentum density [8]. Poynting’s theorem is applied to recover the powerand energy density which are similar if not the same for the majority of for-malisms [10]. Electromagnetic energy density has both microscopic, We,0,and macroscopic, We,1, formsWe,0 =12(0E ·E + 1µ0B ·B) and We,1 = 12(E ·E + µH ·H). (2.7)While energy density and Poynting’s vector, Se, have rarely been disputed,six debated electrodynamic formalisms present unique momentum densityand/or momentum flux density tensors, as shown in Tables 2.4-2.5 [5, 8, 23].In these tables, the variables E = (E · E) 12 , D = (D ·D) 12 , H = (H ·H) 12 ,and B = (B ·B) 12 . The Helmholtz formalism also includes permittivity andpermeability of a Clausius Mossotti medium, such that polynomial expres-sions represent the permittivity,  = r0, and permeability, µ = µrµ0, of amagnetodielectric medium [5].Table 2.4: Electrodynamic Momentum and Power Density.Name Momentum Density Power DensityGe SeMinkowski (MN) [13] D×B E×HAbraham (AB) [14] (E×H)/c2 E×HEinstein & Laub (EL) [15] (E×H)/c2 E×HHelmholtz (HLM) [5] ((µr − 1)(2r + 1) + 3) (E×H)/(3c2) E×HChu (CHU) [19, 20] (E×H)/c2 E×HAmperian (AMP) [19] 0(E×B) E×B/µ0142.2. Electrodynamic sub-system theoryTable 2.5: Electrodynamic Momentum Flux Density Tensors.Name Momentum Flux Density Tensor, TeMinkowski (MN) [13] (D ·E + B ·H)¯I/2−DE−BHAbraham (AB) [14] [(D ·E + B ·H)¯I−DE−BH−ED−HB]/2Einstein & Laub (EL) [15] (0E2 + µ0H2 )¯I/2−DE−BHHelmholtz (HLM) [5](0(2 + 2r − 2r)E2 + µ0(2 + 2µr − µ2r)H2)I/6−EE−µHHChu (CHU) [19, 20] (0E2 + µ0H2 )¯I/2− 0EE− µ0HHAmperian (AMP) [19] (0E2 +B2/µ0 )¯I/2− 0EE−BB/µ0Each of the electrodynamic formalisms can be classified according to theformulation of Maxwell’s equations invoked in their derivation:− Class 1: AB and MN formalisms, constructed from the Minkowskiformulation of Maxwell’s equations (Tab. 2.1).− Class 2: EL and HLM formalisms, constructed from the Einstein &Laub formulation of Maxwell’s equations (Tab. 2.2).− Class 3: CHU and AMP formalisms, constructed from the E-H for-mulation of Maxwell’s equations (Tab. 2.3).Each formulation of Maxwell’s equations is consistent with a Lorentz trans-formation of the reference frame [20], which is a necessary condition forthe kinetic force density, momentum density, and stress predicted by a the-oretical electrodynamics formalism. Recently, both the form of the elec-tromagnetic momentum density and the relativistic frame-invariance of thestress-energy-momentum (SEM) tensor have been actively debated. Shep-pard & Kemp contend that the symmetric Abraham and Einstein & LaubSEM tensors are Lagrangian-variant, while the asymmetric Minkowski SEMtensor is non-causal [64]. Partannen & Tulkki note that the energy con-tinuity equation used by Sheppard & Kemp, which omits the momentumand rest energy of the kinetic sub-system, is incomplete [30]. Present theo-retical debates for a single electrodynamic theory may benefit from generalnumerical methods capable of modeling their observable kinetics [8, 23, 28].152.3. Kinetic sub-system theory2.3 Kinetic sub-system theoryWhile James C. Maxwell is most often credited as the founder of clas-sical electromagnetics, he also founded a dynamical theory of gases [65].Alternative to the phenomenological description of fluids presented by theNavier-Stokes equation, James C. Maxwell’s theory considers an ideal gas asa population of colliding particles that conserve both energy and momentum.Within Maxwell’s description, a probability density function, n(x,v, t), de-scribes the post-collision continuous random velocity, v, of a particle pop-ulation in a multi-dimensional particle velocity space at time t and vectorposition x. In velocity space the zeroth, first, second and third momentsof the probability density function describe the mean mass-density field,ρ(x, t), kinetic momentum-density field, Gk(x, t) = ρ(x, t)u(x, t), kineticenergy density field, Wk(x, t) = ρ(x, t)e(x, t), and the kinetic momentumflux density tensor, Tk(x, t), of an ideal gas, respectively. The vector fieldu(x, t) represents the local mean velocity and the scalar field e(x, t) repre-sents the local energy density per mass density. Each of these moments arelisted in Table 2.6, where∫d3v ≡ ∫ ∫ ∫ dv1dv2dv3.Table 2.6: Moments of the particle population, n(x,v, t).Moment State-Variable Expression at (x, t)0 Mass Density ρ =∫n d3v1 Momentum Density Gk =∫nv d3v2 Energy Density Wk =12∫n|v|2 d3v3 Momentum Flux Density Tensor Tk =∫n(v · v)v d3vMaxwell also defined a probability density function describing the par-ticle population of an ideal gas in a state of equilibrium [65],neq(x, |v˜|, t) = ρ(12piRT)d/2e−|v˜|2/(2RT )∣∣∣∣∣(x,t), (2.8)where v˜(x, t) is the random relative-velocity space, d is the number of di-mensions, R is a gas constant, and T (x, t) is the mean temperature fieldderived from an equation-of-state (EOS). Relative velocity is the differencebetween the random particle population velocity, v(x, t), and the local meanvelocity, u(x, t),v˜(x, t) = v(x, t)− u(x, t). (2.9)162.3. Kinetic sub-system theoryGiven an ideal gas, R = kB/mp, where kB is the Boltzmann constant,and mp is the particle mass. Maxwell’s equilibrium population density,neq(x, |v˜|, t), demonstrates the following properties:1. It is isotropic in velocity space around v(x, t) = u(x, t),2. It is separable as the product of one-dimensional equilibrium distribu-tions [65], and3. It has the same moments of density and energy as n(x, |v˜|, t).Maxwell’s dynamical theory of gases has a philosophical advantage overphenomenological theories of gases; it uses probability axioms to conservemass, energy and momentum while describing their local relative likelihoodat different velocities. Shortly after, Ludwig Boltzmann developed a ki-netic theory of gases by considering the stochastic interactions of ideal-gasparticle populations in velocity space as they approach the equilibrium pop-ulation described by Maxwell [66]. With incorporation of a kinetic force den-sity, fk(x, t), representing an equation of state for a non-ideal gas (liquid),and/or the interaction of external sub-systems (electrodynamic and thermo-dynamic), it has proved useful in the numerical lattice-Boltzmann methodswhich simulate the kinetics of generic fluids [67]. According to Boltzmann’stheory, the kinetic sub-system can be described by Boltzmann’s equation(BE),∂ n∂ t+ v · ∇n+ fkρ· ∂ n∂ v∣∣∣∣∣(x,t)=(dndt)coll∣∣∣∣∣(x,t)= Ω∣∣∣∣∣(x,t). (2.10)Boltzmann’s equation can be derived from the total differential of the parti-cle population density with respect to time, or using the Liouville operator orHamiltonian. The subscript “coll” stands for a collision process as the ordi-nary derivative of the particle population with respect to time, Dt = dn/dt,and is typically described using a collision operator Ω(x, t) [68]. Integra-tion of the collision operator throughout the entire velocity space is zero,∫Ω(x, t)d3v = 0, which conserves both particle populations and their mo-mentum. In general, a collision operator must satisfy the second law of ther-modynamics or Boltzmann’s H-theorem [46]. In Boltzmann’s H-theorem,the quantityH(x, t) =∫n lnnd3v∣∣∣∣∣(x,t), (2.11)172.3. Kinetic sub-system theorycan only decrease, and reaches its minimum value when n(x, t) = neq(x, t)[66]. In Boltzmann’s kinetic theory, as opposed to continuum theory, conser-vation of mass and momentum are related to the collision and movement ofparticle populations via the total probability axiom. In Chu’s E-H formula-tion of electrodynamics and the mass-polariton theory [20, 30], such particlepopulations are used to describe a classical kinetic sub-system. However,in general, Boltzmann’s equation is non-relativistic and its traditional formmay not account for Lorentz transformation or inertial frames of referenceapproaching the speed of light. Development of the relativistic Boltzmann’sequations and lattice-Boltzmann methods for relativistic hydrodynamics arean active area of research outside the scope of this thesis [49, 50].2.3.1 Mass and momentum conservationIn order to demonstrate mass and momentum conservation, vector-spacemoments of a particle-population density are substituted into Boltzmann’sequation to represent particle collisions as a stochastic process. This can beaccomplished by distinguishing zero and non-zero moments of the particle-population density and its derivatives [68]. All moments listed in Table2.6 are non-zero, while any moments of the collision operator, Ω, are zero.Therefore, the equation of mass-conservation,∂∂ t∫nd3v +∂∂ xα∫vαnd3v +Fαρ∫∂ n∂ vαd3v =∫Ωd3v∂ ρ∂ t+∂ ρuα∂ xα= 0,(2.12)can be derived given that∫Ωd3v =∫dndtd3v = 0 and∫∂ n∂ vαd3v = 0. (2.13)Similar derivations can be used for the momentum and energy conservationequations [68],∂ Gk,α∂ t+∂ Tk,αβ∂ xβ= fk,α,∂ Wk∂ t+∂ Sk,α∂ xα= Φk,(2.14)where Gk,α = ρuα are components of momentum density, and Φk = fk,αuα isthe power flux density described using index notation for vector components.182.4. Numerical methodsIf the particle velocity is described using the mean and relative velocity,v = u + v˜, the momentum flux density tensor isTk,αβ = ρuαuβ +∫v˜αv˜βnd3v, (2.15)where the second term describes the stress tensor,σαβ = −∫v˜αv˜βnd3v. (2.16)Since the stress tensor is unknown, the momentum conservation equationis not closed, such that an approximation of n is required by establishinga relation to hydrodynamics via Chapman-Enskog theory (CET). CET isa framework for deriving constitutive relations, such as the Navier-Stokesequation for hydrodynamics, from the molecular parameters described byBoltzmann’s equation. CET establishes a direct relationship between macro-scopic phenomena and microscopic phenomena such as inter-particle colli-sions using an equation-of-state (EOS) [68, 69].2.4 Numerical methodsSeveral numerical methods are used in this thesis including the Yee-lattice finite-difference time-domain method (YL-FDTDM), the electrody-namic lattice-Boltzmann method (ELBM) and the Shan-Chen pseudopo-tential lattice-Boltzmann method (PLBM). Each of these numerical meth-ods discretize the phenomenological equations, or the lattice-Boltzmannequation (LBE), by using central differences or multi-range lattice-stenciloperations. While the discretization of a phenomenological equation isspecific to the represented sub-system, the LBE operations and their dis-cretization are general and can be applied to multiple sub-systems. How-ever, alternative lattice stencils and their weighted operations are requiredto maintain the stability and rotational isotropy of Boltzmann’s equationfor a kinetic sub-system [45]. First this section describes the central dif-ference equation, followed by one-dimensional three-vector (D1Q3), two-dimensional five-vector (D2Q5), two-dimensional nine-vector (D2Q9), three-dimensional seven-vector (D3Q7) and three-dimensional twenty-seven vec-tor (D3Q27) spaces and their general operators. Finally, this section con-cludes with a description of the LBE and its generic algorithm using a two-dimensional kinetic sub-system as an example.192.4. Numerical methods2.4.1 The central difference equationThe central difference equation is a second-order-accurate numerical so-lution for the first derivative, f (1)(x), of a function f(x). central differencescan be derived by subtracting the backward difference Taylor series fromthe forward difference Taylor series of the first derivative of a function. Theforward difference Taylor series isf(x+∆x) = f(x)+∆xf (1)(x)+∆x2f (2)(x)2!+∆x3f (3)(ζ+)3!, ζ+ ∈ (x, x+∆x),while the backward difference Taylor series isf(x−∆x) = f(x)−∆xf (1)(x)+∆x2 f(2)(x)2!−∆x3 f(3)(ζ−)3!, ζ− ∈ (x, x−∆x).The subtraction of these two Taylor series leads to a second-order O(∆x2)central difference approximation of f ′(x),f (1)(x) =f(x+ ∆x)− f(x−∆x)2∆x− ∆x212(f (3)(ζ+) + f(3)(ζ−)). (2.17)The error of the central difference equation is proportional to ∆x2. The one-dimensional central difference equation can approximate partial derivativesin multiple independent/orthogonal dimensions.2.4.2 Lattice stencilsWhile an infinite vector space describes the random relative-velocity inBoltzmann’s equation, only a finite vector set, vm, can exist on a latticedescribed by the LBE [70–72]. Consequently, the lattice-Boltzmann method(LBM) is a finite-difference method in a vector space, as opposed to a scalarspace. Each lattice stencil is described by its number of orthogonal spatialdimensions, and its number of independent discrete lattice-vectors. Thesetwo properties of the lattice stencil determine its order of accuracy, its order-of-isotropy and the phenomena it can describe. In the case of inviscid flow,two-dimensional five lattice vector (D2Q5) and three-dimensional seven lat-tice vector (D3Q7) lattice stencils are sufficient for describing the vectorspace. However, there are general conditions a lattice stencil must satisfyto describe viscous interactions [45, 73]. In addition to mass and momen-tum conservation, rotational isotropy of the lattice is required. In a ki-netic sub-system describing viscous flow, all moments of the lattice-vectorweights, wm ≡ w(|vm|), require isotropy up to the fifth order (as described202.4. Numerical methodsby Maxwell’s theory) [45, 65, 73]. In this thesis, two-dimensional nine latticevector (D2Q9) and three-dimensional twenty-seven lattice vector (D3Q27)stencils are used for their increased stability and accuracy. Both the D2Q9and D3Q27 have the greatest number of lattice vectors that can be repre-sented in a single square lattice unit. Each lattice stencil is illustrated inFig. 2.2.(a) D2Q5 (b) D2Q9(c) D3Q7 (d) D3Q27Figure 2.2: (a) Two-dimension five-velocity (D2Q5) lattice; (b) Two-dimension nine-velocity (D2Q9) lattice; (c) Three-dimension seven-velocity(D3Q7) lattice; (d) Three-dimension twenty-seven velocity (D3Q27) lattice.The black dot represents the present position on the lattice and the nullvector.212.4. Numerical methodsGiven a kinetic sub-system (subscript k), the weight for each discretevelocity in each lattice-stencil, wk,m ≡ wk(|vm|), can be derived from thefollowing conditions [45, 65, 73],Nm∑m=1wk,m = 1,Nm∑m=1wk,mvmα = 0,Nm∑m=1wk,mvmαvmβ = c2sδαβ,Nm∑m=1wk,mvmαvmβvmγ = 0,Nm∑m=1wk,mvmαvmβvmγvmµ = c4s(δαβδγµ + δαγδβµ + δαµδβγ),Nm∑m=1wk,mvmαvmβvmγvmµvmν = 0,(2.18)where cs is the lattice speed of sound representing an isothermal equationof state, and Nm is the number of lattice vectors. Given an electrodynamicsub-system (subscript e), the weight for each non-zero lattice-vector in eachlattice-stencil is we,m ≡ we(|vm| > 0) = 1/4, while the weight for the nulllattice-vector is zero (w(|0|) = 0), which is derived from the constraintsof Eq. (3.37) in Chapter 3 [39]. Table 2.7 describes the lattice stencils interms of their lattice-vectors (with magnitudes |vm|), their lattice-vectorscorresponding weights, and whether they are applied in a kinetic (k) orelectrodynamic (e) sub-system.Table 2.7: Properties of lattice-stencils used in this thesisLattice vm # |vm| wk,m we,mD1Q3[39] 〈0〉 1 0 - 0〈±1〉 2 1 - 1/4D2Q5[39] 〈0, 0〉 1 0 - 0〈±1, 0〉, 〈0,±1〉 4 1 - 1/4D2Q9[74] 〈0, 0〉 1 0 4/9 -〈±1, 0〉, 〈0,±1〉 4 1 1/9 -〈±1,±1〉 4√2 1/36 -D3Q7[39] 〈0, 0〉 1 0 - 0〈±1, 0, 0〉, 〈0,±1, 0〉, 〈0, 0,±1〉 6 1 - 1/4D3Q27[74] 〈0, 0, 0〉 1 0 8/27 -〈±1, 0, 0〉, 〈0,±1, 0〉, 〈0, 0,±1〉 6 1 2/27 -〈±1,±1, 0〉, 〈±1, 0,±1〉, 〈0,±1,±1〉 12√2 1/54 -〈±1,±1,±1〉 8√3 1/216 -222.4. Numerical methods2.4.3 Discrete spatial operators for the lattice-BoltzmannequationDiscrete spatial operators used in the lattice-Boltzmann method can bederived from the lattice stencil, including the velocities and their weightslisted in Table 2.7 [75]. Two common operators used throughout this thesisinclude the gradient and divergence operators,∇ψ(x) =Nm∑m=1w(|vm|)ψ(x + vm)vm=Nm∑m=1wmψ(x + vm)vm,(2.19)∇ · u(x) =Nm∑m=1w(|vm|)u(x + vm) · vm=Nm∑m=1wmu(x + vm) · vm,(2.20)which reduce to the typical one-dimensional central difference equationsgiven a D1Q3, D2Q5, or D3Q7 lattice. Other lattice stencils used, suchas the D2Q9 and D3Q27 lattice stencils, will have more general central dif-ference equations accounting for more directions in the lattice. Multi-rangeoperations may also be used on the lattice, and can be accounted for byincreasing the velocity space to D2Q9r or D3Q27r, where r = 2, 3, ... repre-sents the lattice range [40]. Below is an expansion of Eq. (2.19) for a D2Q9lattice where i and j are the spatial indices in x and y,(∂ ψ∂ x)ij=13(ψi+1j− ψi−1j)+112(ψi+1j+1− ψi−1j+1+ ψi+1j−1− ψi−1j−1),(∂ ψ∂ y)ij=13(ψ ij+1− ψij−1)+112(ψi+1j+1− ψi+1j−1+ ψi−1j+1− ψi−1j−1).(2.21)Given a single lattice and general operators such as Eqs. (2.19) and (2.20)available for lattice-Boltzmann methods, it is rare to find discrete represen-tations such as Eq. (2.21). It is likely such notation is even discouraged,especially since it does not aid in implementing generic lattice-Boltzmannprograms. MATLAB functions that use the notation of Eq. (2.19) for anarbitrary nearest-neighbour lattice-stencil are provided in Appendix D.232.4. Numerical methods2.4.4 The lattice-Boltzmann equationThe lattice-Boltzmann equation (LBE) with no external forces considersthe case of an ideal particle population with a finite number of discretecollision velocities, Nm [76],∂∂ tnm(x, t) + vmα∂∂ xαnm(x, t) = Ω(nm(x, t)), m = 1, ..., Nm, (2.22)where vmα is component α of the lattice vector vm, whilenm(x, t) ≡ n(x,v|v ∈ {vm}, t) (2.23)is a discrete probability mass function equivalent to the conditional probabil-ity that the continuous particle population n(x,v, t) collides and propagateswith a finite number of velocities, v ∈ {vm}, that can be represented on adiscrete lattice. The LBE is generic, such that it can be applied to rep-resent multiple physical sub-systems, such as electrodynamic and kineticsub-systems [37, 39, 76]. However, traditionally it was designed for thekinetic sub-system. There are a number of schemes that have been devel-oped for incorporating an external force density, and different formulationsof the collision operator, which will be described for the kinetic sub-systemin Chapter 4. Moments of nm(x, t) or its equilibrium condition, neqm(x, t),represent different macroscopic state variables in each sub-system. One ex-ample is for a kinetic sub-system, where the zeroth moment is mass density,ρ(x, t), and the first moment is the momentum density, ρ(x, t)u(x, t), whichare obtained using sums as opposed to continuous integrals [71],ρ(x, t) =Nm∑m=1nm(x, t) =Nm∑m=1neqm(x, t),ρ(x, t)u(x, t) =Nm∑m=1nm(x, t)vm =Nm∑m=1neqm(x, t)vm.(2.24)At each discrete time step, nm(x, t) is relaxed as it approaches its equilibriumstate, neqm(x, t), proportional to the numerical order of accuracy. Given akinetic sub-system, orthogonal Hermite polynomials are used to derive andapproximate Maxwell’s equilibrium particle population density to a secondorder of accuracy [76],neqm(x, t) = wmρ(1 +vmαuαc2s+uαuβ(vmαvmβ − c2sδαβ)2c4s)∣∣∣∣∣(x,t), (2.25)242.4. Numerical methodswherecs =1√3δxδt(2.26)is the lattice speed of sound for dimensionless spatial unit δx and dimension-less time unit δt in D2Q9 and D3Q27 lattices. Otherwise, if all higher-orderterms are included in a generic form of the Maxwell equilibrium particle pop-ulation, an entropic relaxation approach is required to accurately representBoltzmann’s H-theorem [46].2.4.5 Collision and propagation operators in thelattice-Boltzmann equationThe lattice-Boltzmann equation is represented in a discrete velocityspace, although its implementation is represented in discrete space and time.Typically, this requires dimensionless spatial unit δx = 1, and a lattice-stencil dependent time unit δt. Particle population nm(x, t) is forced toapproach neqm(x, t) in population collision and propagation operations.The collision operation represents the perturbation of nm(x, t) as itapproaches a state of equilibrium, neqm(x, t), during interparticle collisionsat time t and vector position x. Given a finite number of lattice velocitiesand single relaxation, the discrete Bhatnagar-Gross-Krook (BGK) collisionoperator has been used to represent the relaxation of the particle populationnm(x, t) as it approaches neqm(x, t) according to the discrete Navier-Stokesequation [77],Ωm(x, t) =δtτ(neqm(x, t)− nm(x, t)) , (2.27)where τ is the discrete relaxation time and δt = 1 for a kinetic lattice-Boltzmann method. In practice, the collision operator is added as a differ-ence to nm(x, t) from the previous time step [71], such thatn∗m(x, t) = nm(x, t) + Ωm(x, t),= nm(x, t) +δtτ(neqm(x, t)− nm(x, t)) .(2.28)Following its elastic collision, post-collision particle population n∗m(x, t) isassigned to adjacent lattice node m in the propagation operation.252.4. Numerical methodsThe propagation operation represents the redistribution of n∗m(x, t) toadjacent lattice node m for particle and momentum conservation [76], suchthatnm(x + vmδt, t+ δt) := n∗m(x, t), (2.29)where “:=” is an assignment operator. The collision and consecutive prop-agation process is illustrated in Fig. 2.3. Chapter 4 later discusses entropicand multi-relaxation approaches.(a) Collide (b) Propagate87 8512349( , )mn tx*( , )mn tx ( , )m mn t t t + +x vFigure 2.3: [Color online] Illustration of the collision process on a D2Q9lattice. (a) The particle population outcomes stream toward each latticenode, indexed m = 9, from their respective adjacent lattice nodes, indexedm = 1, ..., 8, and collide elastically. (b) Following the collision, the particlepopulation’s post-collision outcomes propagate to their respective adjacentlattice nodes.26Chapter 3Comparing finite-differencemethods andlattice-Boltzmann methodssolutions of kinetic forcedensityIn considering light as a collisionless gas of pseudo-particles, Maxwell’sequations can be recovered from Boltzmann’s equation [37, 39, 49, 78]. Thischapter demonstrates how such a construction provides time-domain sim-ulations of electric and magnetic fields at the same location and time, in-creasing the numerical accuracy while decreasing the computational time ofan electrodynamic simulation compared with the Yee-lattice finite-differencetime-domain method (YL-FDTDM). First, the initial condition of the elec-trodynamic sub-system, a continuous Hermite-Gaussian wave-function, isdescribed in detail. Second, the YL-FDTDM and discretization of Maxwell’sequations for computing central difference force density are described. Third,the formulation of central difference force-densities is described for Hauser& Verhey’s electrodynamic lattice-Boltzmann method (ELBM) [39]. Thesedescriptions are followed by the implementation methodology of Hauser &Verhey’s ELBM for solving central difference force densities in uniform di-electric media [39], and a comparison between the computational time andnumerical accuracy of force densities evaluated in the YL-FDTDM and theELBM at the same spatial resolution.3.1 Initial conditionOnly monochromatic Hermite-Gaussian beams are used in this thesisbecause: (1) they can be modeled using a paraxial approximation, (2) theyhave finite cross-sectional area, (3) they mimic experimental conditions as-273.1. Initial conditionsociated with focused laser beams. The source condition is constrained to asingle wavenumber and angular time-frequencyk =2piℵλ0and ω =2piT, (3.1)where ℵ is the refractive index of the source medium for free-space wave-length λ0 and time-period T . Free-space wavelength and time-period arelinearly related by the vacuum speed of light, c = 299792458 m/s, such thatλ0 = cT. (3.2)Unless otherwise stated, the source has a refractive index of ℵ = 1 such thatthe wavenumber k = k0 = 2pi/λ0.3.1.1 Hermite Gaussian wave-functionsHermite-Gaussian wave-functions have an amplitude envelope describedby a Gaussian function with a mode order of m along the primary axis, andn along the secondary axis in a transverse plane [79]. In this thesis, the fun-damental TEM00 mode is explored since it represents the output of lasersused in optical radiation pressure experiments which focus the electric fieldto a maximum intensity for maximum optical radiation pressure [22]. TheHermite-Gaussian wave-functions are initialized as periodic and monochro-matic sources, such that they require the minimum number of time samplesfor calculating a time-averaged kinetic force density. The complete Hermite-Gaussian wave-function also has a defined focal position without the use ofoptical lenses [7]. Explicit Hermite-Gaussian solutions can be derived fromthe paraxial Helmholtz equations in Ref. [79], where the wave-functionφmn(r) =√2piw2(z)H(m)(√2xw(z))H(n)(√2yw(z))...× exp[−i(m+ n+ 1)θ(z) + iks22q(z)] (3.3)can represent the field components of both linearly and circularly-polarizedGaussian beams using “physicist’s form” Hermite polynomials H(m) or H(n),which have orders m and n for each transverse dimension. In the generalcase of a three-dimensional Gaussian beam, displacement from the beam-axis along the direction of propagation is represented by z − z0, while s2 =(x− x0)2 + (y − y0)2, i is the imaginary unit, θ(z) = tan−1[(z − z0)/zR)] is283.1. Initial conditionthe Gouy phase, and w(z) = w0√1 + [(z − z0)/zR]2 is the beam-waist fora Rayleigh range zR = kw20/2. The corresponding beam radius, R(z), andcomplex beam parameter, q(z), are defined asR(z) = z − z0 + z2Rz − z0 and1q(z)=1R(z)+ i2kw2(z). (3.4)In three dimensions, displacement is represented by r = (x − x0)xˆ + (y −y0)yˆ + (z − z0)zˆ, where the position of the beam-axis at the beam-waist is(x0, y0, z0). The authors of Ref. [79] used the properties of Hermite polyno-mials to derive Hermite-Gaussian functions for the electric field componentsdescribing linearly and circularly-polarized transverse Hermite-Gaussian, orTEMmn, modes.The linearly-polarized TEMmn Hermite-Gaussian functions for each ofthe electric field components areE(mn)1 (r) = Amnφmn(r),E(mn)2 (r) =Amn4(kw0)2[4mnφm−1,n−1(r)− 2mφm−1,n+1(r)...−2mφm+1,n−1(r) + φm+1,n+1(r)],E(mn)3 (r) =iAmn√2kw0[2mφm−1,n(r)− φm+1,n(r)].(3.5)In three-dimensions the xˆ-polarized TEM00 Gaussian beam can be initializedusing the x-componentE(00)x (s, θ, z) =A00w(z)√2piexp(−iθ(z) + iks22q(z))exp(−iωt). (3.6)However, in two dimensions the direction of propagation is in-plane, suchthat further considerations are necessary. One useful modification is anarbitrary in-plane direction of propagation given by the axial vector s = (x−x0) cosϑxˆ + (y − y0) sinϑyˆ, such that the in-plane curvature of the beam isrepresented by the perpendicular vector s⊥ = (x−x0) sinϑxˆ−(y−y0) cosϑyˆ,where ϑ is the angle of the beam-axis with respect to x. Therefore, in two-dimensional simulations it is appropriate to initialize the Hermite-Gaussianfunction using Hz for a p-polarized simulation asH(00)z (s, θ) =1ηA00w(s)√2piexp(−iθ(s) + iks2⊥2q(s))exp(−iωt), (3.7)293.1. Initial conditionor Ez for an s-polarized simulation,E(00)z (s, θ) =A00w(s)√2piexp(−iθ(s) + iks2⊥2q(s))exp(−iωt), (3.8)where η is the impedance of the source medium. Although a Hermite-Gaussian wave-function of a single component only approximates a solu-tion to Maxwell’s equations, its initialization satisfies the discrete Maxwell’sequations at each consecutive time step to represent a monochromatic sourceof focused electromagnetic radiation. Therefore, simulating beam focusingfrom a lens is not required, such that computational cost and time can bereduced compared with previous studies [7]. Cross-polarization of a linearly-polarized Gaussian beam occurs naturally in three dimensions, such that athree-dimensional simulation of optical radiation pressure from a linearly-polarized Gaussian beam simultaneously represents both predominantly p-polarized and s-polarized electric fields depending on the plane of incidence[79].The circularly-polarized Hermite-Gaussian functions for each of theelectric field components are,E(mn)1 (r) =Amn√2φmn(r),E(mn)2 (r) = iAmn√2φmn(r),E(mn)3 (r) = iAmn2kw0[2mφm−1,n(r)− φm+1,n(r)+ 2inφm,n−1(r)− iφm,n+1(r)],(3.9)for a general TEMmn mode. The circularly-polarized TEM00 mode Hermite-Gaussian beam that propagates in the zˆ direction can only be representedin three dimensions, and requires the initialization of the xˆ and yˆ electricfield components,E(00)x (s, θ, z) =A00√piw(z)exp(−iθ(z) + ik s22q(z))exp (−iωt) ,E(00)y (s, θ, z) =iA00√piw(z)exp(−iθ(z) + ik s22q(z))exp(−iωt).(3.10)This thesis typically defines the field intensity using its peak value, E0, as themaximum electric field solved throughout the spatial domain and temporalduration of the simulation.303.2. The Yee-lattice finite-difference time-domain method3.2 The Yee-lattice finite-difference time-domainmethodThe Yee-lattice finite-difference time-domain method (YL-FDTDM) hasbeen implemented in one, two and three dimensions for solving a wide rangeof electrodynamic phenomena [8, 33, 80]. It is one of the most common andpopular numerical methods in electrodynamics [23, 80]. The YL-FDTDMhas been applied in two dimensions for solving kinetic force densities inponderable media, which has enabled predictions of optical radiation pres-sure [23]. This section describes the construction of central difference forcedensities in the YL-FDTDM, beginning with the evaluation of the electricand magnetic fields, the space and time interpolation, and the momentumdensity and stress tensor [7, 8]. It concludes with a description of numericaldispersion in the YL-FDTDM [80].3.2.1 Electric and magnetic fieldsIn a 2D plane of incidence, there is only an in-plane spatial variationof Maxwell’s equations in Tab. 2.1. Eliminating the third dimension en-ables independent treatment of p- and s-polarized electric fields. The 2Dp-polarized (transverse-magnetic) Faraday and Ampere’s laws are∂ Hz∂ t=− 1µ(σMHz +∂ Ey∂ x− ∂ Ex∂ y),∂ Ex∂ t=1(−σEEx + ∂ Hz∂ y),∂ Ey∂ t=1(−σEEy − ∂ Hz∂ x),(3.11)while the s-polarized (transverse-electric) Faraday and Ampere’s laws are∂ Ez∂ t=1(−σEEz + ∂ Hy∂ x− ∂ Hx∂ y),∂ Hx∂ t=− 1µ(σMHx +∂ Ez∂ y),∂ Hy∂ t=− 1µ(σMHy − ∂ Ez∂ x).(3.12)central differences of these partial-differential equations can be implementedby half-grid space and time offsets of the electric and magnetic field compo-nents. This results in the Yee-lattice structure, which allows numerical so-lutions for electric and magnetic fields in tandem [80]. The two-dimensional313.2. The Yee-lattice finite-difference time-domain methodYee-lattices representing p- or s-polarized planes of incidence are illustratedin Figure 3.1.xyi-1 i i+1j-1jj+1 zHExE yj+12j-12(a) TM, p-polarized Yee-LatticezEHxH yxyi-1 i i+1j-1jj+1j+12j-12(b) TE, s-polarized Yee-LatticeFigure 3.1: In (a), the Yee-lattice representing a p-polarized electric field isillustrated, while in (b), the Yee-lattice representing an s- polarized electricfield is illustrated.The p-polarized partial-differential equations in Eq. (3.11) can be dis-cretized using central differences in space and time to represent the followingcentral difference equationsEn+1xi,j+12=e1xe2xEnxi,j+12+∆t∆yi,j+121e2x(Hn+ 12zi,j+1 −Hn+ 12zi,j),En+1yi+12 ,j=e1ye2yEnyi+12 ,j− ∆t∆xi+12 ,j1e2y(Hn+ 12zi+1,j −Hn+ 12zi,j),Hn+ 12zi,j =h1zh2zHn− 12zi,j −∆tµi,jh2z(Enyi+12 ,j− Enyi− 12,j∆x...−Enxi,j+12− Enxi,j− 12∆y),(3.13)with coefficientse1x = 1−(∆tσEi,j+12/2i,j+ 12), e2x = 1 +(∆tσEi,j+12/2i,j+ 12),e1y = 1−(∆tσEi+12 ,j/2i+ 12,j), e2y = 1 +(∆tσEi+12 ,j/2i+ 12,j),h1z = 1−(∆tσMi,j/2µi,j), h2z = 1 +(∆tσMi,j/2µi,j).(3.14)323.2. The Yee-lattice finite-difference time-domain methodSimilarly, the s-polarized partial differential equations in Eq. (3.12) arediscretized in space and time to represent the following central differenceequations:Hn+1xi,j+12=h1xh2xHnxi,j+12− ∆t∆yµi,j+121h2x(En+ 12zi,j+1 − En+ 12zi,j),Hn+1yi+12 ,j=h1yh2yHnyi+12 ,j+∆t∆xµi+12 ,j1h2y(En+ 12zi+1,j − En+ 12zi,j),En+ 12zi,j =e1ze2zEn− 12zi,j +∆ti,je2z(Hnyi+12 ,j−Hnyi− 12,j∆x...−Hnxi,j+12−Hnxi,j− 12∆y),(3.15)with coefficientsh1x = 1−(∆tσMi,j+12/2µi,j+ 12), h2x = 1 +(∆tσMi,j+12/2µi,j+ 12),h1y = 1−(∆tσMi+12 ,j/2µi+ 12,j), h2y = 1 +(∆tσMi+12 ,j/2µi+ 12,j),e1z = 1−(∆tσEi,j/2i,j), e2z = 1 +(∆tσEi,j/2i,j),(3.16)where ∆x is the dimensional spatial unit and ∆t is the dimensional timeunit.3.2.2 Energy density and momentum densityInterpolation in space and time is required in order to evaluate energydensity and momentum density from the electric and magnetic fields, dueto their offset components in an orthogonal Yee-lattice. The number ofaverages and their construction depends on the number of dimensions in theYee-lattice [7, 8]. In a two-dimensional Yee-lattice representing p-polarizedelectromagnetic fields, interpolation of the electric field components in spaceand time are required to situate these components on the magnetic fieldcomponent to calculate the energy density and momentum density. At timestep n+12 , this is given byEn+ 12xi,j =14(En+1xi,j+12+ En+1xi,j− 12+ Enxi,j+12+ Enxi,j− 12),En+ 12yi,j =14(En+1yi+12 ,j+ En+1yi− 12 ,j+ Enyi+12 ,j+ Enyi− 12 ,j).(3.17)333.2. The Yee-lattice finite-difference time-domain methodThe corresponding energy density isWn+ 12i,j =12(i,jEn+ 12xi,j En+ 12xi,j + i,jEn+ 12yi,j En+ 12yi,j + µi,jHn+ 12zi,j Hn+ 12zi,j), (3.18)while the momentum density is represented asGn+ 12e,xi,j = En+ 12yi,j Hn+ 12zi,j ,Gn+ 12e,yi,j = −En+ 12xi,j Hn+ 12zi,j .(3.19)In a two-dimensional Yee-lattice representing s-polarized electromagneticfields, interpolation of the magnetic field components in space and timeare required to situate these components on the electric field component tocalculate the energy density and momentum density. At time step n+12 , thisis given byHn+ 12xi,j =14(Hn+1xi,j+12+Hn+1xi,j− 12+Hnxi,j+12+Hnxi,j− 12),Hn+ 12yi,j =14(Hn+1yi+12 ,j+Hn+1yi− 12 ,j+Hnyi+12 ,j+Hnyi− 12 ,j).(3.20)The corresponding energy density is,Wn+ 12i,j =12(i,jEn+ 12zi,j En+ 12zi,j + µi,jHn+ 12xi,j Hn+ 12xi,j + µi,jHn+ 12yi,j Hn+ 12yi,j), (3.21)while the momentum density is represented as,Gn+ 12e,xi,j = −En+ 12zi,j Hn+ 12yi,j ,Gn+ 12e,yi,j = En+ 12zi,j Hn+ 12xi,j .(3.22)While these interpolations minimize numerical error in the YL-FDTDM,they still contain numerical error from the dispersion of electric and mag-netic fields in free space. Offset electric and magnetic field components alsorepresent a non-exact magnetodielectric boundary, and introduce numericalreflections. In three-dimensions, they become more involved.3.2.3 Force densityOptical force density in the YL-FDTDM requires intricate spatio-temporalinterpolation and central difference equations that are unique for each elec-tric field polarization. This section is limited to a description of a p-polarized343.2. The Yee-lattice finite-difference time-domain method(TM) electric field Yee-lattice evaluation of force density. The transversemagnetic field solved at time step n+12 is averaged in time and the in-planeelectric-field components at time step n are averaged in space,Hnzi,j =12(Hn+ 12zi,j +Hn− 12zi,j),Enxi,j =12(Enxi,j+12+ Enxi,j− 12),Enyi,j =12(Enyi+12 ,j+ Enyi− 12 ,j).(3.23)The divergence of the momentum flux density tensor,∇ ·Tne,i,j = xˆ∆Tne,xi,j + yˆ∆Tne,yi,j , (3.24)is then solved using a central difference method in space with the interpo-lated field-quantities from Eq. (3.23). Divergence of the Chu momentumflux density tensor in the xˆ-direction is(∇ ·T)ne,xi,j = ... (3.25)− 02∆x((Enxi+1,j)2 − (Enxi−1,j)2) ...+02∆x((Enyi+1,j)2 − (Enyi−1,j)2) ...+µ02∆x((Hnzi+1,j)2 − (Hnzi−1,j)2) ...02∆y(Enxi,j+1Enyi,j+1 − Enxi,j−1Enxi,j−1),while divergence of the Chu momentum flux density tensor in the yˆ-directionis(∇ ·T)ne,yi,j = ... (3.26)02∆x(Enxi+1,jEnyi+1,j − Enxi−1,jEnxi−1,j)...+02∆y((Enxi,j+1)2 − (Enxi,j−1)2) ...− 02∆y((Enyi,j+1)2 − (Enyi,j−1)2) ...+µ02∆y((Hnzi,j+1)2 − (Hnzi,j−1)2) .353.2. The Yee-lattice finite-difference time-domain methodThe derivation for the divergence of each formalisms momentum flux densitytensor is provided in the Appendix B. Appendix D also describes a programimplementation of conventional central difference equations used to evaluatethese force densities. The transfer of kinetic momentum density at time stepn+12 is then calculated via a re-arranged central difference equation of Eq.(2.2),gn+ 12xi,j = gn− 12xi,j − (∆Tne,xi,j )∆t− (Gn+ 12e,xi,j −Gn− 12e,xi,j ),gn+ 12yi,j = gn− 12yi,j − (∆Tne,yi,j )∆t− (Gn+ 12e,yi,j −Gn− 12e,yi,j ),(3.27)followed by the force density at time step n, fnk,yi,j , using a central differencein timefnk,xi,j =1∆t(gn+ 12e,xi,j − gn− 12e,xi,j ),fnk,yi,j =1∆t(gn+ 12e,yi,j − gn− 12e,yi,j ).(3.28)The force density is then averaged over a number of time-periods NT . Givenan angular frequency, ω, the number of computational time steps averagedis Nt = 2piNT /ω∆t and therefore〈fk,xi,j 〉 =1NtNt∑n=0fnk,xi,j ,〈fk,yi,j 〉 =1NtNt∑n=0fnk,yi,j .(3.29)In order to reduce discrete-time error from the irrational time-period, thisthesis recommends NT > 5 [7, 8]. After time-averaging the kinetic forcedensity over NT > 5, there is less than 3% difference in the net force eval-uated. This thesis uses NT = 25 to mitigate time-averaging error in theYL-FDTDM.363.2. The Yee-lattice finite-difference time-domain method3.2.4 Boundary conditionsAbsorbing/radiating boundary conditions are ideal for a time-averagedelectrodynamic simulation. Several different types of these boundary condi-tions exist based on the principle of attenuating free-space matched electricand magnetic field intensity in a small layer of several spatial units. One suchlayer in the YL-FDTDM is the Berenger perfectly matched layer (BPML)[62]. The BPML attenuates both the electric field and the magnetic field ina medium using electric and magnetic conductivity with relative magnitudesequivalent to the free-space impedance, η0,σMµ0=σE0and σM = η20σE . (3.30)The BPML is designed to absorb obliquely incident electromagnetic wavesby attenuating transverse electric and magnetic field components which areoriented tangent to its boundary. Longitudinal field components are notexplicitly attenuated, but indirectly from the transverse field components,which results in numerical reflections from the boundary of the BPML. Nu-merical reflections and dispersion can also occur within the BPML. WhileBerenger had minimized these reflections, a typical BPML requires a thick-ness of at least 3-4 spatial units [62].3.2.5 LimitationsThe YL-FDTDM can solve accurate electric and magnetic fields in thetime-domain in a staggered Yee-lattice provided appropriate space and time-discretization. Since it has well-understood sources of numerical error, theYL-FDTDM can be quickly applied to achieve accurate numerical solu-tions for a wide variety of electromagnetic problems [80]. However, theYL-FDTDM demonstrates several limitations for electrodynamics:− Numerical dispersion in free space, which leads to spurious time-averagedforce densities where there are no electric charges or dipoles,− Space and time offsets between electric and magnetic fields which mustbe interpolated for calculating momentum density and energy density,− Space and time offsets at boundaries, which cause errors in local mo-mentum conservation, leading to numerical reflection and asymmetry,− Indirect relation between electric and magnetic fields to their energyand momentum conservation, which can introduce numerical instabil-ity.373.2. The Yee-lattice finite-difference time-domain methodThe limitations of the YL-FDTDM, and its implementation complexity, haslimited most simulations of kinetic force density from electric and magneticfields to one or two dimensions [7, 8]. However, natural electromagneticradiation is three or four-dimensional (space and time), which has beenconsidered in only a few studies to date [16, 23, 81]. This is largely dueto the implementation complexity and computational time required from aYee-lattice configuration.Numerical dispersion can be described as any deviation of a numericalvalue from its analytical value for a particular frequency bin in a discretefrequency domain. Numerical dispersion occurs naturally in all electromag-netic finite-difference methods as a consequence of finite-difference approx-imations. In the one-dimensional conventional YL-FDTDM, numerical dis-persion is represented bysin(12ω∆t)∆t=1√µ00sin(12k∆x)∆x, (3.31)which can satisfy ω = kc in free space for the condition ∆t = ∆x/c, whichprevents numerical dispersion. In dielectric media or in higher dimensions,however, numerical dispersion cannot be avoided [80]. For example, thetwo-dimensional numerical dispersion relation[1c∆tsin(12ω∆t)]2=[1∆xsin(12kx∆x)]2+[1∆ysin(12ky∆y)]2(3.32)cannot directly satisfy k2x+k2y = ω2c2 in free space [39, 80]. The consequenceof this numerical dispersion in the YL-FDTDM is distortion of the time andspatial periodicity of electromagnetic waves at further distances from thesource, even in the absence of ponderable media. It is necessary that kineticforce density is evaluated everywhere with similar spatial periodicity and ac-curacy, a condition that cannot be satisfied with the frequency-distortion ofelectromagnetic waves. While bandwidth can be increased using a decreasedtime unit [82], spatial unit, higher-order finite-difference [83, 84], or modelorder reduction and eigenvalue perturbation [85, 86], these approaches offera compromise between computational time and accuracy or vice versa. Al-ternatively, the electrodynamic lattice-Boltzmann method (ELBM) avoidsfree-space numerical dispersion in multiple dimensions by including addi-tional degrees of freedom [39]. Fig. 3.2 illustrates the propagation of anelectromagnetic pulse from a point source initialized in both YL-FDTDMand the ELBM. The initial electric field conditionEy(x, y) = 2e−0.05[(x−Nx/10)2+(y−Ny/2)2] (3.33)383.2. The Yee-lattice finite-difference time-domain methodis re-used from Ref. [39] to demonstrate the successful implementation ofthe method and the comparative numerical dispersion of the YL-FDTDM.In the ELBM, δt = 1/2, while in the YL-FDTDM, δt = 1/√2. Eq. 3.33generates an electromagnetic pulse that propagates through a lattice withNx lattice units in the xˆ-direction and Ny lattice units in the yˆ-direction.Fig. 3.2 compares the cross-section of the dimensionless electric field alongits propagation axis, as evaluated in the (a) YL-FDTDM and the (b) ELBM.Spatial distortion in the YL-FDTDM is visible as high-frequency transientsat the end of the pulse which increase over time in Fig. 3.2(a). No such spa-tial distortion is visible in the ELBM Fig. 3.2(b) [39]. The spatial distortionis an effect of numerical dispersion in the YL-FDTDM.(a) YL-FDTDM (b) ELBMFigure 3.2: Comparison of the electric field from an electromagnetic pulsein (a) the YL-FDTDM, and (b) the ELBM during its propagation throughfree space in a dimensionless spatial domain. The plot is a cross-sectionof the two-dimensional simulation along the center of the pulse. Numericaldispersion is evident in (a), yet absent in (b). These sumulations wereconducted by Cael Warner and mimic the simulations described by Hauser& Verhey, 2019 [39].393.3. The electrodynamic lattice-Boltzmann method3.3 The electrodynamic lattice-BoltzmannmethodAlternative to the YL-FDTDM is the recent electrodynamic lattice-Boltzmann method (ELBM) [37, 39, 49, 78]. The ELBM uses a single latticein which electric and magnetic field components all share a common dis-crete position and discrete time. The ELBM derived by Hauser & Verheyuses a single-lattice configuration accomplished by mapping all orthogonalpropagation outcomes of electric and magnetic fields to a single probabilitymass function [39]. The ELBM uses more degrees of freedom to representthe electric and magnetic fields, which provides numerical solutions withgreater accuracy and stability than the YL-FDTDM [39]. Since interparti-cle collisions do not occur, and there is zero viscosity, a simple D3Q7 latticeis sufficient to represent all field propagation outcomes [39]. Collision andpropagation operations in the stochastic Boltzmann process also preservethe space and time period of the electromagnetic waves, with less numericaldispersion which only manifests as a decrease in the field amplitude beyondthe Courant condition. The following sections will describe the advantagesof the ELBM over the YL-FDTDM for electrodynamic simulation.3.3.1 Electric and magnetic fieldsThe ELBM assumes electric and magnetic fields can be represented byan inviscid, collision-less pseudo-particle population density, gm(x, t), as alocal probability mass function which has set of Nm outcomes, where eachoutcome (m = 1, ..., Nm) represents a discrete direction of propagation onthe lattice [37, 39, 49, 78]. Moments of gm(x, t) in terms of orthogonallattice vectors vm, hm, and em represent the power density, magnetic fieldintensity, and electric field intensity in free space [39]. The ELBM operatesat the limit of inviscid flow, τ = 1/2, to represent the linear Ampere andFaraday laws [37, 49]. Given this constraint, the dimensionless time unit isδt = 1/2, and the dimensionless spatial unit is δx = 1. The collision andpropagation operations occur at dimensionless time t and vector position xaccording tog∗m(x, t) = 2geqm (x, t)− gm(x, t)gm(x + vmδt, t+ δt) := g∗m(x, t),(3.34)where “:=” is the assignment operator and g∗m(x, t) is the post-collisionpseudo-particle population [39]. The equilibrium pseudo-particle population403.3. The electrodynamic lattice-Boltzmann methoddensity, geqm (x, t), has been derived using suitable Lyapunov functions toguarantee stability [39, 78]. The equilibrium psuedo-particle populationdensity is represented bygeqm (x, t) =14(E(x, t) · em + H(x, t) · hm), (3.35)where E(x, t) and H(x, t) are the dimensionless electric and magnetic fieldintensity, which are calculated as moments of gm(x, t),E(x, t) =Nm∑m=1gm(x, t)em and H(x, t) =Nm∑m=1gm(x, t)hm. (3.36)The possible lattice vectors depend on the number of dimensions, and in twodimensions, the polarization of the electric field. The lattice vectors can beorganized into matrices, which are referred to as a lattice-matrices in thisthesis. The elements of the lattice-matrix can be derived from the followingconditions for orthogonality in the Chapman-Enskog expansion,Nm∑m=1vαm =Nm∑m=1hαm =Nm∑m=1eαm = 0,Nm∑m=1hαmhβm =Nm∑m=1eαmeβm = 4δαβ,Nm∑m=1eαmhβm = 0,Nm∑m=1vαmeβmhγm = 2αβγ ,Nm∑m=1vαneβmeγm =Nm∑m=1vαmhβmhγm = 0,(3.37)where {α, β, γ} are independent non-equivalent dimensional indices [39].Hauser & Verhey’s method accounts for both left-handed and right-handedorthogonal propagation bases to maximize stability. However, the Minkowskiformulation of Maxwell’s equations in Tab. 2.1 only represents right-handedpropagation bases. Consequently, Hauser & Verhey’s Chapman-Enskog ex-pansion requires a dimensionless time unit of δt = 1/2 [39], to representAmpere and Faraday laws12∇×H = E˙ and 12∇×E = −H˙. (3.38)413.3. The electrodynamic lattice-Boltzmann methodThe factor of 1/2 also manifests in the initial condition for the dimensionlesselectric and magnetic field intensity.3.3.2 Two-dimensional lattice-matricesTwo-dimensional lattice-matrices representNm = 8 non-zero right-handedand left-handed orthogonal bases constructing the outcomes of gm(x, t). Thetwo-dimensional lattice-matrices for a p-polarized, or transverse-magnetic(TM), electric field are[vαm] =[1 0 -1 0 -1 0 1 00 -1 0 1 0 1 0 -1],[eαm] =[0 -1 1 0 0 -1 1 0-1 0 0 1 -1 0 0 1],[hαm] =[-1 -1 0 0 1 1 0 0].(3.39)Fig. 3.3 illustrates the right-handed bases which represent half the propaga-tion outcomes of a p-polarized electric field. The remaining left-handed baseshave an opposite order and orientation. The polarization density, P(x, t),and magnetization density, M(x, t), are stationary outcomes of E(x, t) andH(x, t).12341234(a) (b) (c)mv me mh1234P MFigure 3.3: [Color online] Illustration of discrete lattice vectors formingright-handed propagation bases for a p-polarized electric field on a D2Q5 lat-tice. The remaining lattice vectors form the left-handed propagation baseswhich have an opposite order and orientation on the lattice. Polarizationand magnetization density represent stationary electric and magnetic fieldoutcomes. The symbol “⊗” is “into-the-page,” while the symbol “ ” is“out-of-the-page.”423.3. The electrodynamic lattice-Boltzmann methodThe two-dimensional lattice-matrices for an s-polarized, or transverse-electric (TE), electric field are[vαm] =[1 0 -1 0 -1 0 1 00 -1 0 1 0 1 0 -1],[eαm] =[1 -1 0 0 1 -1 0 0],[hαm] =[0 1 -1 0 0 -1 1 0-1 0 0 1 1 0 0 -1].(3.40)Fig. 3.4 illustrates the right-handed bases as half the propagation outcomesof an s-polarized electric field. The remaining left-handed bases have anopposite order and orientation. The polarization density, P(x, t), and mag-netization density, M(x, t), are stationary outcomes of E(x, t) and H(x, t).(b)1234(a) (c)mv me mh12341234PMFigure 3.4: [Color online] Illustration of discrete lattice vectors that formthe right-handed propagation bases for an s-polarized electric field on aD2Q5 lattice. The remaining lattice vectors form the left-handed propaga-tion bases which have an opposite order and orientation on the lattice. Polar-ization and magnetization density represent stationary electric and magneticfield outcomes. The symbol “⊗” is “into-the-page,” while the symbol “ ”is “out-of-the-page.”433.3. The electrodynamic lattice-Boltzmann method3.3.3 Three-dimensional lattice-matricesThe three-dimensional lattice-matrices represent Nm = 12 non-zero left-handed and right-handed orthogonal bases constructing the outcomes ofgm(x, t),[vαm] =1 0 0 -1 0 0 -1 0 0 1 0 00 -1 0 0 1 0 0 1 0 0 -1 00 0 -1 0 0 1 0 0 1 0 0 -1 ,[eαm] = 0 0 -1 0 1 0 0 0 -1 0 1 0-1 0 0 0 0 1 -1 0 0 0 0 10 -1 0 -1 0 0 0 -1 0 -1 0 0 ,[hαm] = 0 1 0 0 0 -1 0 -1 0 0 0 10 0 1 -1 0 0 0 0 -1 1 0 0-1 0 0 0 -1 0 1 0 0 0 1 0 .(3.41)These lattice-matrices are later used in a three-dimensional multi-physicssimulation of optical radiation pressure. Two-dimensional and three-dimensionalimplementations of the ELBM are the same, except for their lattice-matrices.Therefore, a numerical method solving a three-dimensional force density issimpler to implement in the ELBM than in the YL-FDTDM.3.3.4 Ponderable mediaIn ponderable media, the propagation of electric and magnetic fields areimpeded by electric and magnetic dipoles at time t and vector position x.The corresponding dimensionless polarization density, P(x, t), and magneti-zation density, M(x, t), manifest as the outcomes for electric and magneticfields remaining at their present locations [78]. These outcomes are forcedto approach their equilibrium conditions, Peq(x, t) = (r(x)− 1)E(x, t) andMeq(x, t) = (µr(x)− 1)H(x, t), with a collision operation only (no propaga-tion since they are stationary)P(x, t+ δt) = 2Peq(x, t)−P(x, t),M(x, t+ δt) = 2Meq(x, t)−M(x, t). (3.42)In the presence of ponderable media, the electric and magnetic fields musthave re-normalized moments of gm(x, t) to reflect the outcomes of polariza-tion or magnetization density [37, 39, 78]. The dimensionless electric and443.3. The electrodynamic lattice-Boltzmann methodmagnetic fields in lossless magnetodielectric media are thereforeE(x, t) =1r(x)Nm∑m=1(gm(x, t)em + P(x, t)),H(x, t) =1µr(x)Nm∑m=1(gm(x)hm + M(x, t)).(3.43)While Eq. (3.43) provides a natural and accurate description of a ponderablemedium’s magnetization and polarization density, the description requiresadditional operations to represent a point charge within a single lattice unit[39].3.3.5 Energy density and momentum densityThe energy density is the same for each of the electrodynamic formalisms,and can therefore be represented as the centred moment of geq,2m (x, t),We(x, t) =Nm∑m=1geq,2m (x, t). (3.44)When the dimensionless velocity is chosen as c0 = 1, both the dimensionlesspower density and the momentum density are the first moment of geq,2m (x, t)in terms of vm,Ge(x, t) =Nm∑m=1geq,2m (x, t)vm. (3.45)This suggests an energy continuity equation of the form 12∇·Se = −W˙e [39],which is true for uniform non-dispersive magnetodielectric media. Such anenergy continuity condition requires that components of E(x, t) or H(x, t)are initialized at half their desired value to remain wholly consistent withMaxwell’s equations. Provided such a scaling in the initial condition, the di-mensional conversions areW′e(x′, t′) = We(x′, t′)0, S′e(x′, t′) = Se(x′, t′)/η0,and G′e(x′, t′) = Ge(x′, t′)/(c2η0) remain consistent with the ELBM Faradayand Ampere laws in Eq. (3.38) and Tab. 2.1. Each of the electrodynamicdimensional conversions are summarized in Tab. 3.1.453.3. The electrodynamic lattice-Boltzmann method3.3.6 Force densityThe force density is the sum of time variation in momentum density andspatial variation in the momentum flux density tensor. The momentum fluxdensity tensor divergence can be represented using a single-lattice divergenceoperator that is the same independent of the number of lattice dimensions.In the case of the D1Q3, D2Q5 and D3Q7 lattices used in Hauser & Ver-hey’s ELBM, the divergence of Chu’s momentum flux density tensor can berepresented as(∇ ·Te(x, t)) = ...181δxNm∑m=1[(|E(x + vˆmδx, t)|2 + |H(x + vˆmδx, t)|2)I...− 2 (E2(x + vˆmδx, t) + H2(x + vˆmδx, t)) ]vˆm,(3.46)where vˆm represents the unit vector of lattice velocity vm, and the dimen-sional scaling is(∇ ·T′e(x′, t′)) =0∆x(∇ ·Te(x′ , t′)) (3.47)since Eq. (3.46) reduces to a central difference equation in a D2Q5 or D3Q7lattice. Meanwhile, the dimensional time-variation of energy density andmomentum density can be respectively obtained using central differenceequations∂∂ tW′e(x′, t′) =121δt∆t(W′e(x′, t′+ δt∆t)−We(x′ , t′ − δt∆t))and∂∂ tG′e(x′, t′) =121δt∆t(G′e(x′, t′+ δt∆t)−Ge(x′ , t′ − δt∆t)).(3.48)When calculated at each full time step, or every two iterations in the two-dimensional or three-dimensional ELBM, Eqs. (3.47) and (3.48) retain theirsecond order of accuracy. The dimensional kinetic force density can then becalculated as the negative sum of Eqs. (3.47) and (3.48)f′k(x′, t′) = −∇ ·T′e(x′, t′)− ∂∂ tG′e(x′, t′). (3.49)Since only a single lattice is used in the ELBM, space and time interpolationare no longer necessary to evaluate this kinetic force density.463.3. The electrodynamic lattice-Boltzmann method3.3.7 Dimensional conversionsTab. 3.1 provides the dimensional conversions for the ELBM. Thesedimensional conversions were derived from the equations provided by Hauser& Verhey in their method, and in the programmatic implementation of theirmethod [39]. The dimensional conversions listed in Tab. 3.1 require that anelectric or magnetic field component is initialized at one-half of its desiredamplitude. Tab. 3.1 is a subset of the larger table of conversions providedin Tab. C.1 of Appendix C.Table 3.1: Conversion of dimensionless electrodynamic sub-system variablesVariable Name Dimensionless (x, t) Dimensional (x′, t′)Time t t′= t∆t/2Position vector x x′= x∆xElectric field intensity E E′= EMagnetic field intensity H H′= H/η0Energy density We W′e = We0Poynting’s vector Se S′e = Se/η0Momentum density Ge G′e = Ge/(c2η0)Momentum-flux tensor Te T′e = Te03.3.8 LimitationsUnlike in the YL-FDTDM, the ELBM time- and spatial-units cannot bemodified with respect to one another since the relaxation time is restricted toa single value, τ = 1/2. The consequence is that adaptive mesh refinementmaybe impractical. Therefore, there maybe limited systems that benefitfrom the ELBM’s increased accuracy given its fixed space-time relation ina uniform lattice. However, static or sub-optimal dynamic mesh refinementin the YL-FDTDM can lead to late-time instability from increasing numer-ical error at interfaces between different spatial grids [87]. Therefore, theconstraint of a uniform mesh with square lattice units is often used in elec-trodynamic simulations focused on predicting optical radiation pressure [7].473.4. Comparative simulation configuration3.4 Comparative simulation configurationBoth the spatial unit, ∆x, and the time unit, ∆t, are simultaneouslyand proportionally varied to represent the same spatio-temporal domain.The source is time harmonic with a wavelength of λ = 530 nm and a time-period of T = 1.77 fs in free space. The instantaneous kinetic force densitiesare averaged over NT = 25 optical periods. Two slightly different systemsare required to compare the two methods. Fig 3.5 illustrates the 5λ × 5λYL-FDTDM spatial domain, and the 4λ×4λ ELBM spatial domain for eval-uating force densities. The two domains cannot be exactly identical, sincethe YL-FDTDM requires additional spatial padding (shaded gray region)for its combined BPML (λ/4 thickness) and reflection buffer (λ/4 thick-ness). The reflection buffer is necessary to avoid time-averaged numericalreflections in free space from the source and BPML. However, in both illus-trated systems, the black region is exactly the same size (4λ× 4λ) with thesame spatial resolution. It is within this region that the computational timeand accuracy of the YL-FDTDM and ELBM time-averaged kinetic forcedensities are compared.Although the YL-FDTDM requires less computational time than theELBM for evaluating electric and magnetic fields at the same spatial reso-lution, the YL-FDTDM requires greater computational time for evaluatingelectric and magnetic fields of similar accuracy. Hauser & Verhey noted thatthe YL-FDTDM’s computational time is 18 times greater than the ELBM’scomputational time in achieving solutions of the same numerical accuracy[39]. This thesis compares the additional computational time required forthe YL-FDTDM and ELBM to solve time-averaged force density in the blackregion, which is the same spatial domain with the same spatial resolutionfor the two configurations illustrated in Fig. 3.5. Similarly, the numeri-cal accuracy of the net force in the black region is compared between theYL-FDTDM and ELBM.483.4. Comparative simulation configuration44(b)D1.333=5445(a)D 1.333=Figure 3.5: Illustration of the spatial domain used in (a) the YL-FDTDM,and in (b) the ELBM. Each spatial domain has an equivalent number oflattice units per wavelength (resolution), although a different size. Half-wavelength, λ/2, boundaries (gray regions) are included as a buffer in theYL-FDTDM for the quarter-wavelength λ/4 Berenger perfectly matchedlayer and its numerical reflections. The source condition is represented asa white dashed line near the bottom of each spatial domain. The blackregion represents where the force density is evaluated; this truncated spatialdomain has the same size in the YL-FDTDM and the ELBM, with the samespatial resolution. The diameter of the dielectric cylinder cross-section isD = 2λ, and its refractive index is ℵ = 1.333 (where r = 1.778 and µr = 1).493.5. Comparative analysis3.5 Comparative analysisThe Hermite-Gaussian beam is initialized using a magnetic field polar-ized in the zˆ direction and an electric field polarized in the xˆ direction,such that electromagnetic waves propagate in the yˆ direction on a two-dimensional xy-plane. Therefore, the electric field is p-polarized and parallelto the plane of incidence. The cylinder cross-section has a relative permit-tivity of r = 1.778 and permeability of µr = 1, and is non-conductive. Theaccuracy of the YL-FDTDM and the ELBM are compared using a time-averaged net line-force density. The net time-averaged line-force density intwo-dimensions can be calculated from〈F〉 =∫A〈f〉dA′ ≈∑i,j〈fi,j〉∆x∆y ≡∑〈x〉=〈i,j〉〈f(x, t)〉∆x∆y, (3.50)which is equivalent for each electrodynamic formalism, given horizontal andvertical spatial indices i and j. In the ELBM, x = 〈i, j〉, since δx = 1 andδy = 1, such that the same numerical integration method is used to evaluatenet time-averaged line-force density in the YL-FDTDM and the ELBM. Dueto symmetry, the analytical solution for the xˆ-directed net time-averagedline-force density is 〈Fx〉 = 0. The xˆ-directed net time-averaged line-forcedensity is used to compare the YL-FDTDM and ELBM for its simplicity.The YL-FDTDM and the ELBM are compared in terms of computationaltime and absolute error with increasing spatial resolution from nλ = 10 tonλ = 80 spatial units per wavelength.3.6 Computational time and accuracycomparisonBoth the YL-FDTDM and the ELBM are used to predict central differ-ence force densities according to multiple electrodynamic formalisms. Simi-lar computational time is required by each electrodynamic formalism, withmedium-independent formalisms requiring the least computational time. Atthe same spatial resolution, the YL-FDTDM is more efficient, such that itcan compute a 50% greater sized array for the electric and magnetic fieldsin the same computational time. However, additional computational time isneeded to evaluate the central difference time-averaged kinetic force density.Fig. 3.8(a) shows that, within the subset domain (black region), the space-and time-interpolated force densities in the YL-FDTDM require 4 to 8 timesgreater computational time than the single-lattice force densities evaluated503.6. Computational time and accuracy comparisonin the ELBM. Fig. 3.8(b) shows that the YL-FDTDM net time-averagedline-force density contains 11 to 35 times greater error than the ELBM’s atthe same spatial resolution. This corroborates an earlier study, where it wasshown that the ELBM evaluated more accurate electric and magnetic fieldsthan the YL-FDTDM at the same spatial resolution [39].(a) YL-FDTDM C.T. (b) YL-FDTDM C.T. of 〈f〉(c) ELBM C.T. (d) ELBM C.T. of 〈f〉Figure 3.6: [Color online] Computational time (C.T.) of the YL-FDTDM(a,b), and the ELBM (c,d) as the number of spatial units per wavelength, nλ,increases. Surplus computational time required to evaluate central differenceforce densities are illustrated in (b) for the YL-FDTDM, and (d) for theELBM. Simulations were performed in MATLAB using an Intel(R) Xeon(R)CPU E3-1225 v3 @ 3.20GHz processor having 7.80 GB RAM.513.6. Computational time and accuracy comparison(a) (b)Figure 3.7: The net time-averaged line-force density in the x-direction (→),〈Fx〉, plotted (a) in natural units and (b) on a semi-logarithmic scale. Ana-lytically, 〈Fx〉 = 0 pN/m.(a) (b)Figure 3.8: (a) Computational time ratio and (b) error ratio plotted withrespect to the number of spatial units per wavelength, nλ. When comparedwith the ELBM, the YL-FDTDM requires 4 to 8 times greater computa-tional time for recovering solutions with 11 to 35 times greater error.523.6. Computational time and accuracy comparison3.6.1 Time-averaged force densityChu’s formalism predicts time-averaged force densities throughout thecross-section of the dielectric cylinder. Comparing Chu’s time-averaged forcedensities in the ELBM and YL-FDTDM at nλ = 10 and nλ = 40, the YL-FDTDM predicts much greater non-physical time-averaged force densitiesin free space (which should not exist in the absence of electric dipoles).(a) YL-FDTDM, nλ = 10 (b) ELBM, nλ = 10(c) YL-FDTDM, nλ = 40 (d) ELBM, nλ = 40Figure 3.9: [Color online] The Chu time-averaged force density in the x-direction (→) for a cylindrical cross-section in (a) the YL-FDTDM and(b) the ELBM, at nλ = 10. Time-averaged force densities in free spaceare greater in the YL-FDTDM than in the ELBM. The YL-FDTDM time-averaged force densities in (c) approach the ELBM time-averaged force den-sities in (d) at nλ = 40.53Chapter 4Multi-physics simulation ofoptical radiation pressure ona multi-phase fluidSince an analysis of deformation is required to discriminate between elec-trodynamic formalisms, this chapter describes an algorithm for the multi-phase multi-physics simulation of optical radiation pressure. In this chapter,we combine the ELBM from Chapter 3 with a lattice-Boltzmann methoddescribing a multi-phase fluid. First this chapter provides a brief literaturereview of several different multi-physics simulations of optical radiation pres-sure on fluids and the advantages and disadvantages of the lattice-Boltzmannmethod. Following this literature review is a description of the differentcomponents of the kinetic lattice-Boltzmann method used in this thesis andhow they comprise the fluid-solver algorithm. The algorithm requires afluid relaxation time that is significantly greater than the time required forelectromagnetic waves to propagate through the fluid medium. This as-sumption is typically true for optical radiation pressure on quasi-stationarymedia, where it has been successfully applied in previous electrodynamicsimulations [7, 8, 23].4.1 Review of multi-physics simulation methodsThe multi-physics simulation of optical radiation pressure on fluids hasbeen accomplished using an analytical approach with a boundary-integral-element method [26, 27], the finite-difference-time-domain (FDTD) methodwith front-tracking [7, 80, 88], and the finite-element method (FEM) [23,28]. Each approach is philosophically similar in that it discretizes the phe-nomenological Navier-Stokes and Maxwell’s equations in space and/or timeto numerically solve the deformation of a ponderable medium. Quantitativecomparisons between experimentally measured and predicted deformationof a fluid-fluid interface are typically sufficient [27, 28]. However, these544.1. Review of multi-physics simulation methodsprevious studies analyzed a solution for the trivial Laplace-pressure, whichmay not distinguish between predictions of optical radiation pressure forseveral different formalisms [7, 8, 29]. Finite-element and front-trackingmethods can also enforce the definition of an interface, which preventsa description of natural multiphase fluid separation, vaporization or coa-lescence under the influence of optical radiation pressure. Finite-elementand front-tracking methods also solve the nonlinear advection term in theNavier-Stokes equation (NSE) non-locally, which requires implicit iterativesolutions for hyperbolic differential equations [43]. Alternatively, the lattice-Boltzmann method (LBM) indirectly solves NSE via Boltzmann’s equationand solves the non-linear advection term locally. However, the greatestadvantages of the LBM are its simplicity and achievement of direct massand momentum conservation via a total probability axiom. Recently, en-tropic relaxation LBMs also have been developed for representing turbulentflow with increased stability [46], multi-range LBMs have been developedfor compressible gases [89], and relativistic LBMs have been developed fordescribing gases with particle velocities approaching the speed of light [58].Most multi-phase fluid simulations using an LBM invoke either a free-energy or interaction-potential (pseudopotential) approach for describingthe phase separation. The pseudopotential approach was developed first[67], but its description of an effective mass density was found thermody-namically inconsistent [40]. Consequently, a thermodynamically consistentfree-energy approach was developed [90]. However, the free-energy approachrequires more computational operations, and its initial implementation ex-hibited Galilean invariance and instability for greater liquid-vapour densityratios [40]. While both the pseudopotential and free-energy approach havebeen modified to overcome their limitations, the pseudopotential approachis related to microscopic interactions while remaining simpler to implement.In its simplest form, the pseudopotential approach is also Galilean invari-ant and can approximate thermodynamic consistency provided a physicalequation of state [40, 91], which is sufficient for a number of fluid experi-ments [40]. Thermodynamic inconsistency of the pseudopotential methodis a consequence of second-order treatment of the third-order term in theShan-Chen force density, the third-order term describing the phase separa-tion [92]. Increased academic interest in the pseudopotential approach hasprovided greater information on adaptations overcoming its thermodynamicinconsistency [40, 92–94], and spurious velocities at an interface [40, 92].Recently, the pseudopotential lattice Boltzmann method (PLBM) was ap-plied to analyze the deformation of multi-phase fluids from an electrostaticfield or continuous optical radiation pressure [36, 38]. Each study assumed554.2. The pseudopotential lattice-Boltzmann methoda single electrodynamic formalism describing electrohydrodynamic forces oroptical radiation pressure at the interface [36, 38]. This thesis applies a sim-ilar single-component multiphase PLBM to simulate the kinetic sub-system[36, 38, 40, 67], but uses the ELBM to model a generalized electrodynamicsub-system [39]. The generalized sub-system interaction enables studies ofoptical radiation pressure from focused Hermite-Gaussian beams represent-ing experimental lasers [79], and their scattering from a multiphase lineardielectric fluid in both two and three dimensions. This thesis also accountsfor the existence of multiple electrodynamic formalisms that may predictunique deformation of fluids in the rest frame [7, 16, 23].4.2 The pseudopotential lattice-BoltzmannmethodThe pseudopotential lattice-Boltzmann-method (PLBM) is a numericalsolver of Boltzmann’s equation for multi-phase microscopic and macroscopic(mesoscopic) non-ideal fluids, which uses a non-ideal force density and equa-tion of state to represent a vapour-liquid interface. The central parameterof the PLBM is a probability mass function (PMF), nm(x, t), representing aparticle population in a random discrete-velocity space, vm, at dimensionlesstime t and vector position x, where m = 1, ..., Nm as an index for the discretevelocities. The PMF, nm(x, t), is made to approach the Maxwell equilibriumprobability mass function (EPMF), neqm(x, t), using either a single-relaxation,multi-relaxation, or entropic-relaxation process. Moments of the PMF in ve-locity space only represent the mass and momentum density of an ideal gas,and require the addition of an external force to represent their non-idealstate.4.2.1 Initial conditionThe initial condition for the multi-phase studies in this thesis is a dis-tributed mass-density in the shape of a superellipse, a cylindrical cross-section, or a sphere with negligible velocity. Regardless of the initial den-sity, the most stable interface is represented by a hyperbolic tangent func-tion which serves as an analytical solution to Maxwell’s construction for anisothermal liquid vapor phase coexistence [95],ρinit.(r, t) =12(ρliq. + ρvap.)± 12(ρliq. − ρvap.) tanh(2(r− ri) · nˆiι), (4.1)564.2. The pseudopotential lattice-Boltzmann methodwhere ri is the locus of the interface with respect to the origin, nˆi is thenormal unit vector, ι is its thickness, and r is the displacement vector froma common origin.The superellipse is initialized as an approximate rectangle,ρinit.(y, t) =12(ρliq. + ρvap.)± 12(ρliq. − ρvap.) tanh(2(y − h)ι), (4.2)where y is the dimension in which mass density varies, and h representsthe liquid depth. When initialized within one half of a macroscopic fluiddomain, the superellipse can stabilize to represent an approximately flat in-terface. Flat interfaces are important since they have been used in formerexperiments and simulations of optical radiation pressure [7, 28]. However,a flat interface is thermodynamically unstable at the microscale accordingto Laplace’s law [96]. Even at millimetric scales, which exceed 1000 opticalwavelengths, the steady-state water interface has a curved meniscus. Non-equilibrium oscillations from the unstable initial condition of a flat interfacecan contaminate the numerical solution. At the very least, additional com-putational time or specialized boundary conditions are required to stabilizea superellipse.The cylinder cross-section or sphere is initialized asρinit.(r, t) =12(ρliq. + ρvap.)± 12(ρliq. − ρvap.) tanh(2(r − ri)ι), (4.3)where r is the radial distance from the center of the droplet, ri is the radiusof the interface representing a cylinder or sphere, and ι is the thickness of theinterface. Adding the second term represents a droplet, while subtractingthe second term represents a bubble. The cylindrical cross-section or spher-ical droplet with negligible velocity has the minimum free surface energyand the greatest stability provided an ambient hydrostatic pressure. There-fore, it represents ideal microscopic light-matter interactions free from anyexternal media. Droplets can also mimic cells which are analyzed in opticaltweezing experiments [52]. Despite the challenges of accurately detectingsuch a system in experiment, it is convenient for testing an electrodynamicsimulation method and its application.574.2. The pseudopotential lattice-Boltzmann method4.2.2 Mass and momentumMass density, ρ(x, t), and momentum density, ρ(x, t)u(x, t), of an idealgas can be computed as moments of the PMF, nm(x, t), in its discretevelocity space [71],ρ(x, t) =Nm∑m=1nm(x, t), ρ(x, t)u(x, t) =Nm∑m=1nm(x, t)vm. (4.4)The momentum density of a non-ideal gas (or liquid), however, is perturbedby an external force density (such as the pseudopotential force density),which must be incorporated [67]. This thesis uses the Shan-Chen forcingscheme [40],u(x, t) =1ρ(x, t)(Nm∑m=1nm(x, t)vm +τδt2f(x, t)), (4.5)and the general Guo forcing scheme [97],u(x, t) =1ρ(x, t)(Nm∑m=1nm(x, t)vm +δt2f(x, t)), (4.6)to represent phase separation and external forces. Guo’s forcing scheme isused in single-relaxation and multi-relaxation approaches for two-dimensionalfluid systems. Meanwhile, Shan-Chen’s forcing scheme is used in single-relaxation approaches for three-dimensional fluid systems.4.2.3 The pseudopotential force densityThe pseudopotential, or Shan-Chen, force density represents a third-order instability between a liquid and its vapor phase [40, 92]. While theinstability has been criticized for its thermodynamic inconsistency, whichmay not allow explicit definition of surface tension, an equation of statemay be used to indirectly approximate thermodynamic consistency and cor-responding surface tension [40, 91]. In this case, the Shan-Chen force densitybased on an effective density is given byfSC(x, t) = −Gψ(x, t)Nm∑m=1w(|vm|) [ψ(x + vmδt, t)] , (4.7)584.2. The pseudopotential lattice-Boltzmann methodwhile the corresponding effective density approximating a thermodynami-cally consistent interface may be represented asψ(x, t) =√2(p(x, t)− c2s)Gc2s, (4.8)where w(0) = 4/9, w(1) = 1/9, and w(2) = 1/36 and the coefficient G = 1.Given the effective mass, ψ(x, t), any arbitrary equation of state can be usedto define p(x, t) [91]. The Carnahan-Starling equation of state demonstratesthe greatest thermodynamic consistency and stability [91],p(x, t) = ρRT 1 + bρ/4 + (bρ/4)2 − (bρ/4)3(1− bρ/4)3 − aρ2∣∣∣∣∣(x,t)(4.9)where a = 1, b = 4, and R = 1 are dimensionless versions of fundamen-tal constants. However, it is also possible to generate a customized cubicequation of state [98]. Consequently, the PLBM has found wide applica-tion and adaptation for representing fluid phenomena including microscopicdroplets and bubbles [40]. The Shan-Chen pseudopotential force density isalso amenable to single-relaxation [40], multi-relaxation [93], and entropic-relaxation LBMs [99]. The Shan-Chen force density is essential to the pseu-dopotential approach and it has been applied in Guo and Lycett-Brown& Luo forcing schemes [92, 100] (among a variety of alternative forcingschemes) which perturb the PMF during the collision operation [40].594.2. The pseudopotential lattice-Boltzmann method(a) (b)Figure 4.1: [Color online] (a) The Carnahan-Starling equation of state (CS-EOS) liquid-vapor phase coexistence curve in terms of reduced temperature(T = T ′/T ′c ) and reduced density (ρ = ρ′/ρ′c). The theoretical predictionis derived by equating chemical potentials in the liquid and vapor phases,similar to Ref. [91]. (b) Dimensionless “[ ]” mass density, ρ, plotted withits maximum liquid-vapour density ratio, ρliq./ρvap. = 1180 at T ′ = 0.531,on a two-dimensional 128×128 grid.4.2.4 The kinetic force densityMultiple force densities can be added to the Shan-Chen force densityvia the property of superposition, such that force densities can be incorpo-rated between interacting sub-systems provided they have equivalent spatio-temporal dimensions. Since each system property is dimensionless, it is pos-sible to linearly scale the multi-physics system to any physical dimensionsprovided that the mean-free-path of the fluid is smaller than the spatial unit.Regardless of the scale, dimensional conversion is required to represent phys-ical system properties. In each sub-system, the dimensionless spatial unitis the fundamental unit of measure, δx = 1, by which mass density andtime, ρ and δt = 1, are secondary units of measure that can be scaled usingthe dimensional two-phase critical point mass density, ρ′c, and physical timeunit, ∆t. The time-averaged kinetic force density is a component of externalforce density added to the Shan-Chen force density and incorporated withinany of the forcing schemes described in this thesis.604.2. The pseudopotential lattice-Boltzmann method4.2.5 Collision operation and operatorsCollisions in the PLBM are locally non-linear, and bypass the Poissonequation with pseudo-compressibility, which enables their simple paralleliza-tion and adaptation for different scales and experiments [41, 42]. Given thatthe collision outcomes are linearly independent, the PLBM operations canbe parallelized for greater computational efficiency in solving unsteady orincompressible fluid flows [42]. However, the original PLBM exhibits thegreatest spurious liquid velocities at an interface [42]. Modifications to thecollision operator include a third-order forcing term with an accurate equa-tion of state [91, 92], multi-relaxation [40, 93], and/or entropic relaxation[46, 99].The single relaxation method uses the Bhatnagar-Gross-Krook (BGK)collision operator [77], with the PMF, nm(x, t), approaching the EPMF,neqm(x, t), according ton∗m(x, t) = nm(x, t)−δtτ(nm(x, t)− neqm(x, t)) ...+(1− δt2τ)Sm(x, t)δt.(4.10)The single-relaxation EPMF can be represented in the time-domain us-ing Chapman Enskog expansion with Hermite polynomials [76]. The termSm(x, t) is a non-ideal forcing term that modifies the collision operator foran external force density according to the forcing scheme. Multiple forc-ing schemes exist, but this thesis applies the Shan-Chen and Guo forcingschemes for their stability or accuracy [40, 97]. However, the Shan-Chen andGuo forcing schemes are limited compared with alternative forcing schemes,such as Lycett-Brown & Luo’s forcing scheme, which has tunable surfacetension and interface width for increased stability and versatility [92]. Inthe Shan-Chen forcing scheme [40],Sm(x, t) = 0, (4.11)which causes the density ratio between the liquid and vapor phases to changewith relaxation time, τ , causing thermodynamic inconsistency and instabil-ity with deviation from τ = 1 [93, 97]. This density-dependent relaxationtime in the Shan-Chen forcing scheme severely limits the PLBM’s range ofviscosity [93]. In this thesis, the Shan-Chen forcing scheme is only usedin qualitative three-dimensional simulations for its stability given a relax-ation time of τ = 1. Alternatively, Guo’s forcing scheme can be used to614.2. The pseudopotential lattice-Boltzmann methodincorporate the Shan-Chen force density [97],Sm(x, t) =(1− 12τ)wm[vm − u(x, t)c2s+vm · u(x, t)c4svm]· f(x, t), (4.12)and offers enhanced thermodynamic consistency with a greater range of re-laxation time or viscosity. While the Shan-Chen forcing scheme compromisesthermodynamic consistency for stability, Guo’s forcing scheme compromisesstability for thermodynamic consistency. Regardless, both forcing schemesare only conditionally stable since they approximate the third-order Shan-Chen force density as second-order [92]. To represent the complete third-order Shan-Chen force density with tunable surface tension, the Lycett-Brown-Luo third-order forcing scheme is recommended [92].The multi-relaxation method has a collision operation in moment spacethat is transformed to velocity space prior to the streaming/propagationoperation [101]. The collision operation is the same in velocity space andmoment spacen∗m(x, t) = nm(x, t)−Nq∑q=1Λmq(nq(x, t)− neqq (x, t))...+Nq∑q=1(δmq − 12Λmq)Sq(x, t)δt,(4.13)nˆ∗m(x, t) = nˆm(x, t)−Nq∑q=1Λˆmq(nˆq(x, t)− nˆeqq (x, t))...+Nq∑q=1(δmq − 12Λˆmq)Sˆq(x, t)δt,(4.14)where δmq is the Kronecker-delta symbol representing elements of the N×Nidentity matrix for a DdQN lattice [93, 101]. This thesis implements theMRT collision operator on a D2Q9 lattice only [93, 102], although alter-native implementations exist in three dimensions for D3Q15 and D3Q19lattices [103, 104]. Multi-relaxation allows independent tuning of relaxationtimes for each element of nˆeqm(x, t) based on bulk and kinematic viscosity,subject to constraints of hydrodynamics and stability [105]. The relaxationparameters are organized using a collision matrix in either velocity space,624.2. The pseudopotential lattice-Boltzmann method[Λ], or moment-space, [Λˆ]. However, the collision matrix in moment spacecan be defined as a diagonal (diag) matrix [105][Λˆ] = diag[s1, s2, s3, s4, s5, s6, s7, s8, s9], (4.15)where the relaxation parameters, s1, ..., s9, modify each component of [nˆ], amatrix defining nˆm(x, t) [93]. When defined in a moment space, the collisionmatrix de-couples the physical quantities of nˆeqm(x, t). Each parameter froms1 to s9 can be tuned independently, subject to stability constraints [105].One such constraint is that s1 = s4 = s6 = 1. The last two parameters, s8and s9, are related to the single relaxation time and dimensionless kinematicviscosity vias8 = s9 =1τand ν =(τ − 12)c2sδt. (4.16)One additional parameter, s2, is related to bulk viscosity, ζ,s2 =(ζc2sδt+12)-1, (4.17)while the remaining parameters s3 and s5 = s7 are independently tuned forstability [93]. The condition s5 = s7 is required to satisfy symmetry [93, 101].Conversions between velocity space and moment space are accomplishedusing transformation matrix, [M]. Given a D2Q9 lattice, velocity-space tomoment-space and moment-space to velocity-space transformations are[nˆ] = [M][n] and [n] = [M]-1[nˆ], (4.18)given transformation matrix [106][M] =1 1 1 1 1 1 1 1 1-4 -1 -1 -1 -1 2 2 2 24 -2 -2 -2 -2 1 1 1 10 1 0 -1 0 1 -1 -1 10 -2 0 2 0 1 -1 -1 10 0 1 0 -1 1 1 -1 -10 0 -2 0 2 1 1 -1 -10 1 -1 1 -1 0 0 0 00 0 0 0 0 1 -1 1 -1. (4.19)The vector representing the moment-space PMF, [nˆ], is[nˆ] =[ρ,Wk,W2k , Gk,x, Sk,x, Gk,y, Sk,y, σk,xx, σk,xy]ᵀ, (4.20)634.2. The pseudopotential lattice-Boltzmann methodwhere ρ is the mass density, Wk is the kinetic energy density, Gk,x and Gk,yare the xˆ and yˆ components of the kinetic momentum density, Sk,x and Sk,yare similar components of the kinetic energy flux density, and σk,xx and σk,xyare the diagonal and off-diagonal components of the stress-tensor represent-ing the kinetic sub-system [93]. In Eq. 4.20, the symbol “ᵀ” representsthe transpose operation. In moment space, as opposed to velocity space,each of the outcomes have a physical meaning. Similarly, the equilibriumdistribution and forcing term can be transformed to moment space using[nˆeq] = [M][neq] and [Sˆ] = [M][S], (4.21)which results in[nˆeq] =ρ−2ρ+G2k/9ρρ−G2k/9ρGk,x−Gk,xGk,y−Gk,y(G2k,x −G2k,y)/ρGk,xGk,y/ρand [Sˆ] =06(uxfk,x + uyfk,y)−6(uxfk,x + uyfk,y)fk,x−fk,xfk,y−fk,y2(uxfk,x − uyfk,y)uxfk,y + uyfk,x, (4.22)when expanded for a D2Q9 lattice [93]. While independent tuning of eachrelaxation time allows the MRT approach to preserve a mass-density ratioover a greater range of kinematic viscosity than the SRT approach, the MRTforcing term requires minor adjustment for thermodynamic consistency andstability with increasing density ratio. Minor adjustment of [Sˆ] in Eq. (4.22)can tune the third-order instability of the Shan-Chen force density for in-creased thermodynamic consistency in the MRT approach [102],[Sˆ] =06(uxfk,x + uyfk,y) +12σγs2ψ2δt(1−0.5s2) |fk|2−6(uxfk,x + uyfk,y)− 12σγs3ψ2δt(1−0.5s3) |fk|2fk,x−fk,xfk,y−fk,y2(uxfk,x − uyfk,y)uxfk,y + uyfk,x, (4.23)644.2. The pseudopotential lattice-Boltzmann methodwhere s2 and s3 are relaxation parameters from Eq. 4.15, ψ is the equiv-alent mass density from Eq. 4.8, |fk| = (f2k,x + f2k,y)12 is the magnitude ofthe kinetic force density, and σγ is a tuning parameter for the third-orderinstability [102]. Such adjustment is necessary to achieve high density ratiostypical of the Carnahan-Starling equation of state. Since this thesis does notrepresent particular fluid experiments, it implements existing MRT relax-ation parameters to demonstrate generic fluid behaviour [93, 102]. Furthermodifications to [Sˆ] are necessary to achieve an adjustable surface tension[103, 104, 107].The entropic-relaxation method has a Karlin-Boltzmann-Collision (KBC)operator that represents the Boltzmann-H theorem [46]. The KBC collisionoperator is distinct from the BGK collision operator in that it maximizesentropy to represent a physical process. Recently, the KBC operator hasbeen applied to free-energy and pseudopotential LBMs for multiphase flu-ids, with numerical solutions of increased accuracy compared with the BGKcollision operator [99, 108]. Since negligible turbulence is demonstrated inthe studies of this thesis, as per their momentum conservation, the KBCoperator was found unnecessary. Nonetheless, this thesis recommends theKBC operator for future studies examining more general fluid behaviour[109].4.2.6 Equilibrium probability mass functionThe Maxwell equilibrium probability mass function (EPMF) is a discreterepresentation of the continuous Maxwell-Boltzmann distribution represent-ing a population density of colliding ideal-gas particles [65]. EPMFs can beapproximately derived for SRT or MRT methods using Hermite polynomials,or exactly derived for an entropic relaxation method using a separation ofthe EPMF in each orthogonal dimension [46, 76]. While each EPMF can beused to represent the equilibrium of the kinetic sub-system investigated inthis thesis, entropic relaxation approaches were only recently developed forthe PLBM and are still being refined for general implementation [99, 109].Therefore, the kinetic second-order discrete Hermite-polynomial EPMF andthe BGK collision operation were found to be sufficient as a proof of concept.654.2. The pseudopotential lattice-Boltzmann methodBhatnagar-Gross-Krook collision operators, either employing sin-gle or multi-relaxation, approximate the EPMF in a velocity space. TheHermite-polynomial EPMF was found to be stable with appreciable accu-racy in both the two-dimensional and three-dimensional multiphase LBMsused in this thesis [76],neqm(x, t) = w(|vm|)ρ(x, t){1 +1c2svm · u(x, t)...+12c4s(vm · u(x, t))2 − 12c2su2(x, t)}.(4.24)However, KBC operators need an exact definition of the EPMF to represententropic relaxation in three dimensions, which has increased stability forrespresenting a wider range of viscosity in optical radiation pressure experi-ments [46, 109]. Such methods maybe useful for investigating photothermalbubble nucleation or vaporization [54].Karlin-Boltzmann collision operators have increased stability and ac-curacy by directly obeying the Boltzmann H-theorem. BGK SRT and MRTcollision operators ignore entropy, as do classical methods representing theNSE, which assume that entropy is a higher-order numerical error. How-ever, this assumption can cause both classical solvers of the NSEs and LBMsto become unstable or lose accuracy as time progresses. Karlin-Boltzmanncollision (KBC) operators maximize entropy to minimize growth in numer-ical error, helping close the LBE for hydrodynamics and combat a loss innumerical accuracy with time [46]. Implementation of the KBC operatoris beyond the scope of this thesis, although it has recently been appliedto model multiphase liquid-vapor water using either a free-energy LBM orPLBM [99, 108].4.2.7 Propagation operationFollowing the collision operation in either the SRT or MRT approach,the post-collision PMF is re-distributed to adjacent lattice nodes, such thatnm(x + vmδt, t+ δt) := n∗m(x, t), (4.25)where “:=” is an assignment operator. This propagation or “streaming”operation is similar to all LBMs.664.2. The pseudopotential lattice-Boltzmann method4.2.8 Boundary conditionsBoundary conditions in kinetic LBMs redistribute outcomes nm(x, t) ata boundary to represent a particular macroscopic fluid behaviour at thatboundary. However, there exist more outcomes nm(x, t) than macroscopicmoments such that an inverse solution (nm(x, t) as a function of macroscopicvariables such as ρ(x, t) or Gk(x, t)) is not unique. Therefore, hundreds ofboundary conditions exist for kinetic lattice-Boltzmann methods as of 2021.This thesis implements periodic, bounce-back, and non-equilibrium bounce-back boundary conditions. Each of these boundary conditions prevent parti-cles from leaving the simulation domain. While periodic boundaries conservemomentum at each instant in time, the bounce-back boundaries do not.Periodic boundaries allow the transfer of mass and momentum out ofone boundary, where they simultaneously re-enter the system through theopposite boundary. Consequently, the entire simulation domain represents asingle spatio-temporal period of an infinite periodic system. Such boundaryconditions are advantageous for simulations of optical radiation pressure,since they conserve mass and momentum at all times. Periodic boundariescan be implemented in lattice Boltzmann methods by assigning the PMFon one boundary to its opposite boundary [110],nm(x, t) := nm(x + L, t), (4.26)which is illustrated in Fig. 4.2. Virtual nodes on each boundary repre-sent the boundary nodes on the opposite side of the simulation domain.Only outcomes of nm(x, t) directed into the virtual space are assigned tonm(x + L, t). In this thesis, periodic boundary conditions are used for sim-ulating the deformation of droplets.581( , )mn tx ( , )mn t+x L581(a)673( , )mn tx ( , )mn t+x L673(b)Figure 4.2: Illustration of periodic boundary conditions674.2. The pseudopotential lattice-Boltzmann methodBounce-back boundaries reflect particles as an opposite case to periodicboundaries, and are expressed asnm(x,vm, t) = nm(x,−vm, t+ δt), (4.27)along the boundary [110]. Fig. 4.3 illustrates the perfect reflection at theboundary, which does not move. Since the bounce-back boundary does notmove following a collision, it does not necessarily conserve momentum [110].Further, the bounce-back boundary is second order accurate with respect tospace [O(δx2)] which can affect the solution of the third-order Shan-Chenforce density and its instability everywhere. Additional momentum exchangealgorithms have been developed to account for the transfer of momentumfrom the fluid to the boundary [111], however, these algorithms can becumbersome compared with a higher-order bounce-back boundary scheme.( )m t t+v( )m tvFluidWallSolidFigure 4.3: Particle with mass m and velocity v reflects from a rigid wall,such that its velocity is reversed. Since v(t+δt) = −v(t), the time-averagedparticle momentum is 〈mv〉 = (mv(t) + mv(t + δt))/2 = 0. However,momentum transferred to the wall becomes 2mv, such that momentum isnot conserved at each moment in time.684.2. The pseudopotential lattice-Boltzmann methodNon-equilibrium bounce-back boundary conditions Alternatively,the non-equilibrium bounce-back boundary is a third order accurate [O(δx3)]boundary scheme for straight walls [112]. It uses the mass density and ve-locity of the fluid at the wall in order to enforce a no-slip condition. If an ad-ditional external force density is acting on the wall, it must also be included[113]. In a two-dimensional simulation using a D2Q9 lattice with veloci-ties index m = 1, ..., 9, each with magnitude |vm| = 1, the non-equilibriumbounce-back boundary conditions have a unique assignment at each bound-ary− West boundary:ρw =11− uw,x[n9 + n2 + n4 + 2(n3 + n6 + n7)− fw,x2],n1 = n3 +23ρwuw,x − fw,x6,n5 = n7 − 12(n2 − n4) + 12ρwuw,y +16ρwuw,x − fw,x6− fw,y4,n8 = n6 +12(n2 − n4)− 12ρwuw,y +16ρwuw,x − fw,x6+fw,y4.(4.28)− East boundary:ρw =11 + uw,x[n9 + n2 + n4 + 2(n1 + n5 + n8) +fw,x2],n3 = n1 − 23ρwuw,x − fw,x6,n7 = n5 +12(n2 − n4)− 12ρwuw,y − 16ρwuw,x +fw,x6+fw,y4,n6 = n8 − 12(n2 − n4) + 12ρwuw,y − 16ρwuw,x +fw,x6− fw,y4.(4.29)− North boundary:ρw =11 + uw,y[n9 + n1 + n3 + 2(n2 + n6 + n5) +fw,y2],n4 = n2 − 23ρwuw,y +fw,y6,n7 = n5 +12(n1 − n3)− 12ρwuw,x − 16ρwuw,y +fw,x4+fw,y6,n8 = n6 − 12(n1 − n3) + 12ρwuw,x − 16ρwuw,y − fw,x4+fw,y6.(4.30)694.2. The pseudopotential lattice-Boltzmann method− South boundary:ρw =11− uw,x[n9 + n1 + n3 + 2(n4 + n7 + n8)− fw,x2],n2 = n4 +23ρwuw,y − fw,y2,n5 = n7 − 12(n1 − n3) + 12ρwuw,x +16ρwuw,y − fw,x4− fw,y6,n6 = n8 +12(n1 − n3)− 12ρwuw,x +16ρwuw,y +fw,x4− fw,y6.(4.31)In equations Eqs. (4.28)-(4.31), uw,x and uw,y represent the components ofthe local mean fluid velocity at the wall/boundary, ρw represents the massdensity at the wall/boundary, and fw,x and fw,y represent the componentsof the force density at the wall/boundary.( )m t t+v( )m tvFluidWallSolid( )tF( )w tuFigure 4.4: [Color online] The particle population with mass m and velocityv reflects from the rigid wall, but its reflection is affected by non-equilibriumvelocity and forces acting at the wall.4.2.9 Dimensionless numbersDimensionless numbers have unique definitions in kinetic LBMs, such asthe PLBM, compared with their definition for measured dimensional values.This is because the lattice speed of sound in two and three dimensions, cs =(1/√3)δx/δt is distinct from the dimensional speed of sound. Consequently,the corresponding dimensionless numbers relevant to this thesis includeRe =|u|nxδxc2s(τ − δt2 )and Ma =|u|cs, (4.32)which are the Reynolds number and Mach number, respectively. The simu-lations in this thesis distinguish between cases of incompressibility Ma < 0.1and weak compressibility Ma ∈ (0.1, 0.3).704.2. The pseudopotential lattice-Boltzmann method4.2.10 Dimensional conversionTable 4.1 provides a list of kinetic dimensional conversions used in thischapter. The dimensional conversions listed in Tab. 3.1 are a subset of thelarger table of conversions provided in Tab. C.1 of Appendix C.Table 4.1: Conversion of dimensionless kinetic sub-system variablesVariable Dimensionless (x, t) Dimensional (x′, t′)Position vector x x′= x∆xTime t t′= t∆tVelocity u u′= u(∆x/∆t)Mass density ρ ρ′= ρρ′cForce density fk f′k = fk(ρ′c∆x/(∆t)2)Momentum density Gk G′k = Gk(ρ′c∆x/∆t)Similar conversions exist for kinematic viscosity, ν, or bulk viscosity, ζ,which can be derived from standard units. Given a BGK collision operation,the dimensional relationship between the spatial unit, ∆x, and time unit,∆t, is often quadratic. The kinematic viscosity (and fluid velocity) for aparticular ∆t and ∆x is decided by relaxation time or parameters. Therelaxation time is related to both the time unit and the kinematic viscosityof the fluid for a single-relaxation method,ν = c2s(τ − δt2). (4.33)Typically, the kinematic viscosity and velocity of the fluid are known, whilethe relaxation time is either τ = 1 or tuned to a desired value to minimizesecond-order error. Therefore, a dimensional form of Eq. 4.27 is used toderive the dimensional time unit,∆t =(∆x)2c2sν′(τ − δt2), (4.34)where cs = 1/√3 in D2Q9 and D3Q27 lattice stencils. Using an example,if the dimensional spatial unit ∆x = 26.5 nm, and the kinematic viscosityof water is ν′= 1.0034 × 10-6 m2/s, the corresponding time unit is ap-proximately ∆t = 1.05 ns. The bulk viscosity can be similarly defined torepresent a physical fluid with increased stability for the multi-relaxationapproach described in section 4.2.5.714.3. Algorithm4.3 AlgorithmFig. 4.5 shows a schematic of the combined electrodynamic and ki-netic lattice-Boltzmann method. While both the electrodynamic and ki-netic sub-systems have an identical spatial unit, ∆x, the electrodynamicsub-system has a negligible time unit compared with the kinetic sub-system,∆te ∆tk, allowing a time-averaged kinetic force density to be computed.The time-averaged kinetic force density is imported into the kinetic sub-system, provided its dimensionless conversion. Perturbation of the fluid isallowed to occur until an appreciable deformation of the median fluid densityby ∆x, which is checked by the “shape change” condition in the kinetic sub-system. If “shape change” is “True”, then the electrodynamic sub-systemis re-initialized, and the time-averaged force density is recomputed and re-applied to the kinetic sub-system to complete one algorithm cycle. Thealgorithm simulates the continuous deformation of a fluid subjected to asteady and continuous time-harmonic Gaussian beam. Time-averaged forcedensities existing in free space, which exist as an error of a finite time aver-age in either the YL-FDTDM or ELBM, are set to zero to prevent numericalinstability in the vapour phase.Mass density and velocityCompute equilibriumPlotIncrement timeShape changePeriodic or bounce-back boundaryPropagateCollideFieldsCompute equilibriumIncrement timeFree-boundaryPropagateCollideInitialize mass densityInitialize electric fieldTime-averaged force density(a) Electrodynamic sub-system (b) Kinetic sub-system,E H 𝑔𝑒𝑞 𝑔 = 2𝑔𝑒𝑞 − 𝑔et tt = + kt tt = +, u eqninit. n n= +FalseTrueFigure 4.5: Multi-physics simulation algorithm for the combined interactionof an (a) electrodynamic sub-system, and (b) kinetic sub-system. In betweenthe two sub-systems is the transfer of a time-averaged kinetic force densityfrom the electrodynamic sub-system to the kinetic sub-system, given that∆te ∆tk.72Chapter 5Two- and three-dimensionaloptical acceleration anddeformation of multi-phasefluidsNumerical simulations have demonstrated that multiple electrodynamicformalisms predict identical rigid-body kinematics [7, 8]. However, such con-clusions have not held for fluidic or viscoelastic media [7, 8, 23, 28, 114]. Onesignificant difference established by the YL-FDTDM with front-tracking isthat the Chu and Amperian/Lorentz formalisms predict concave fluid-fluidinterface deformation from a p-polarized beam [7, 8]. Concave interfacedeformation was predicted by the Chu and Amperian/Lorentz formalismsregardless of surface tension. For example, similar interface deformationswere predicted for illuminated air-water or toluene-water interfaces, eventhough their surface tensions differ by four orders of magnitude [7, 8]. Pre-vious studies which investigated multiple electrodynamic formalisms pre-dictions of optical radiation pressure on macroscopic fluids, which used ei-ther finite-element or front-tracking methods, omitted a study of the fluid’smomentum conservation [7, 28]. These studies were further restricted toconditions of axial symmetry in three dimensions [28], or the condition oflow density ratio in two dimensions [7]. This chapter uses an alternativeapproach, the lattice-Boltzmann method, to represent the deformation ofgeneric multiphase fluid-fluid interfaces subjected to optical radiation pres-sure in both two and three dimensions [36, 38]. Specifically, this chapterinvestigates the interaction between the electrodynamic lattice-Boltzmannmethod (ELBM) described in Chapter 3 and the psuedopotential lattice-Boltzmann method (PLBM) described in Chapter 4 according to the al-gorithm illustrated in Fig. 4.5. Multiphysics simulations are performedin two-dimensional and three-dimensional spatial configurations representa-tive of the optical acceleration and trapping of microscopic quasi-stationary735.1. Two-dimensional simulationmulti-phase fluid interfaces and droplets. This chapter demonstrates that,in support of theoretical arguments [115], there is identical interface defor-mation of dielectric fluids predicted by multiple electrodynamic formalismsin the rest frame. However, these formalisms fall within separate classes asdiscussed in Chapter 2, each class having a unique formulation of Maxwell’sequations [5, 7, 16, 20, 27, 115, 116]. If the classes are treated equivalentlyusing the same multi-physics continuity equations [115], each class may pre-dict distinct deformation of a dielectric fluid.5.1 Two-dimensional simulationIn this section, two-dimensional deformation and optical accelerationfrom focused s- and p-polarized Hermite-Gaussian beams is investigated forthe cases of both a flat liquid-vapour interface and a circular liquid-vapourinterface.5.1.1 MethodologyThe multi-phase fluid is representative of water at 70 ◦C, with a liquidphase of dimensional mass density ρ′liq. = 977 kg/m3 and kinematic viscosityν′= 1.0034 × 10−6 m2/s, but with a different surface tension which can-not be explicitly represented in the Guo or Shan & Chen forcing schemes.Each continuous-wave time-harmonic Gaussian beam is initialized with anoptical wavelength of λ0 = 530 nm, a vertical (yˆ) direction of propagation,and is focused to a beam-waist of w0 = λ0 at the center of the simulationdomain. Refer to Chapter 3 “Initial condition” for a complete descriptionof the Hermite-Gaussian beam profile. Two-dimensional simulation is suffi-cient for analyzing deformation along a single plane of incidence from eitheran s- or p-polarized Gaussian beam that propagates in the yˆ direction withits beam axis intersecting (x0, y0). While the choice of the electrodynamiclattice-Boltzmann method is invariant of the quasi-stationary multi-phasefluid investigated, the pseudopotential lattice-Boltzmann method must beadapted to satisfy constraints of stability, thermodynamic consistency, andnumerical accuracy of a particular fluid system. Therefore, a different adap-tation of the PLBM is described for each fluid simulation.745.1. Two-dimensional simulationThe flat liquid-vapour interface, which is modeled using a single relax-ation (SRT) PLBM with a multi-range lattice (D2Q92), is used to examinethe incompressible and weakly compressible deformation of dielectric liquidwith an optical index of water (ℵ = 1.333) predicted by multiple electrody-namic formalisms. The flat liquid-vapour interface is initialized at y = 0 nmin a spatial domain with a horizontal dimension x ∈ (−2120, 2120) nm anda vertical dimension y ∈ (−1060, 1060) nm. The dimensionless Shan-Chensurface tension of the liquid-vapour interface is measured as γSC = 0.235with determination coefficient R2 = 0.9647. Non-equilibrium bounce-back(NEBB) boundary conditions are initialized along each of the walls (West,North, East, and South) in the spatial domain. Fig. 5.1 illustrates thesimulation configuration representing a two-dimensional flat liquid-vapourinterface.0 0w =North NEBBWestNEBBSouth NEBBykxEastNEBBy3liq.6 21.7791.333977 kg/m' 1.0034 10 m /s'r −==== 3vap.6 2' 0.828s1.001.00kg/m' 1.0034 10 m /r −==== m(2120,1 60) n0m(2120, 1060) n−m( 2120, 1 n060)− −m( 2120,1060) n−(0,0)Figure 5.1: [Color online] The configuration used in a two-dimensional sim-ulation of optical radiation pressure on a flat interface. The dimensions ofthe simulation domain are 8λ0 × 4λ0, bordered by a black line representingthe non-equilibrium bounce-back (NEBB) boundaries. Black arrows fromthe boundary represent the reflection of particle populations. The thin grayline represents the interface, while the gray region represents the dielectricfluid. The white region is free space. The green dashed line is the sourcecondition, which is a Hermite-Gaussian beam with a curvature representedby the thin dashed black lines, and a beam waist of w0 = λ0 (not to scale)centred at the position (0,0).755.1. Two-dimensional simulationThe circular liquid-vapour interface, which is modeled using a multirelaxation (MRT) PLBM with a nearest-neighbour lattice (D2Q9), is usedto examine weakly compressible optical deformation of a dielectric or mag-netic droplet with an optical index of water (ℵ = 1.333) predicted by mul-tiple electrodynamic formalisms. Unlike the SRT PLBM, the MRT PLBMcan be tuned for stability and thermodynamic consistency to model eitherdroplets or bubbles of the same diameter with a liquid-vapor density ratioexceeding 1000 (ρliq./ρvap. > 1000) using the stability adjustment parameterσγ = 0.195. The circular liquid-vapour interface representing a “droplet”with a diameter of, D = 1.5λ0, (λ0 = 530 nm) and refractive index ℵ = 1.333is centred in the xy-plane (x0, y0) = (0, 0). Fig. 5.2 illustrates the simu-lation configuration representing a two-dimensional circular liquid-vapourinterface.North Periodic0 0w =WestPeriodicSouth PeriodicykxEastPeriodicy3liq.6 21.779 or 1.7791.333977 kg/m' 1.0034 10 m /s'r r −= ==== 3vap.6 2' 0s1.001.001.00kg/m' 1.0034 10 /.828mrr −===== m(1060,1 60) n0m(1060, 1060) n−m( 1060, 1 n060)− −m( 1060,1060) n−(0,0)01.5D =Figure 5.2: [Color online] The configuration used in a two-dimensional simu-lation of optical radiation pressure on a circular interface (cylindrical cross-section). The size of the spatial domain is 4λ0 × 4λ0 and bordered by ablack line representing periodic boundaries. Black arrows at the boundariesrepresent propagation of particle populations through the boundaries. Thegray region represents the dielectric or magnetic fluid which has a diame-ter D = 1.5λ0, and properties listed in the bottom right-hand corner. Thewhite region is free space, with properties listed in the top right-hand cor-ner. The green dashed line represents the position of the source, which isa Hermite-Gaussian beam with a curvature represented by the thin dashedblack lines and a beam waist of w0 = λ0 (not illustrated to scale).765.1. Two-dimensional simulation5.1.2 Flat liquid-vapor interfaceIncompressible deformation and flow occurs for a Mach number lessthan 0.1 (Ma < 0.1), as a figure of merit in kinetic lattice-Boltzmann simu-lations. Fig. 5.3(a) shows the deformation of an incompressible (Ma < 0.1)dielectric liquid (r = 1.779) being subjected to optical radiation pressurefrom a p-polarized Hermite-Gaussian beam. The simulations were run upto 370 ns to represent the transient state of the liquid prior to its steadystate, after which the liquid-vapour interface stabilizes within the Hermite-Gaussian beam and no longer exhibits steady deformation (but small peri-odic oscillations). The Abraham (AB), Minkowski (MN), Einstein & Laub(EL), and Helmholtz (HLM) electrodynamic formalisms predict an identi-cal convex (outward) liquid-vapour interface deformation representative ofa positive pressure. The Chu and Amperian/Lorentz (AMP/CHU) electro-dynamic formalisms are identical in a dielectric, and predict the oppositecase: a concave (downward) liquid-vapour interface deformation represen-tative of a negative pressure. Chu and Amperian/Lorentz predictions of aconcave interface deformation are consistent with Refs. [7, 8], which arethe only distinction between electrodynamic formalisms predictions of time-harmonic time-averaged optical radiation pressure on a dielectric in twodimensions. Fig. 5.3(b) shows the maximum deformation over time, whichis distinct between the cases of convex and concave deformation.0 0 0 0Figure 5.3: [Color online] Dielectric liquid-vapor (r = ℵ2 and µr = 1)interface deformation simulated using PLBM given p-polarized incident ra-diation with peak intensity E0 = 0.875 GV/m. The CHU/AMP formalismsare identical and opposite of the remaining formalisms, which are degener-ate.775.1. Two-dimensional simulationWeakly compressible deformation and flow occurs for a Mach num-ber in the range of Ma ∈ (0.1, 0.3). Given weak compressibility, variationsin mass density cause a medium to become inhomogeneous. In this study,the difference in mass density of the liquid phase from the vapor phaseis linearly mapped to the difference in the relative permittivity of waterfrom free space. Therefore, spatial variations in mass density of the liquidphase cause proportional spatial variations in dielectric permittivity. How-ever, mass density and dielectric permittivity are uncorrelated at opticalfrequencies [117]. Their relation is described as electrostriction, which is atheoretical contribution to the deformation of the fluid nanoseconds afterits illumination [27, 28]. Otherwise, electrostriction has balanced surfaceand volume contributions at optical frequencies [118]. Spatial variations inmass density and dielectric permittivity are similar to the spatial variation ofthe bulk force density, which predicts increased fluid deformation comparedwith a surface force density. Therefore, assuming dielectric permittivity islinearly proportional to mass density may lead to discrepancies in fluid de-formation predicted by different electrodynamic formalisms. However, socan the assumption that dielectric permittivity remains homogeneous withvariation in mass density. Although the source of differences in deformationpredicted by bulk force density and surface force density requires furtherinvestigation, they may restrict the accuracy of the simulation method tothe incompressible Mach flow regime (Ma < 0.1).0 0 0 0Figure 5.4: [Color online] Dielectric liquid-vapor (r = ℵ2 and µr = 1)interface deformation simulated using the PLBM, given s-polarized incidentradiation with peak intensity E0 = 1.125 GV/m. Weak compressibilityleads to differences in deformation predicted by the EL/HLM/CHU/AMPformalisms and by the MN/AB formalisms.785.1. Two-dimensional simulation5.1.3 Circular liquid-vapour interfaceFig. 5.5(a) shows the deformation of a two-dimensional weakly compress-ible (Ma < 0.3) dielectric droplet (r = 1.779), represented as a cylindricalcross section (or circular interface), subjected to a p-polarized Gaussianbeam with a peak electric field intensity of E0 = 0.98 GV/m. The AB, MN,EL, and HLM electrodynamic formalisms predict similar elongation of thedroplet longitudinal to the direction of propagation. CHU and AMP for-malisms predict identical elongation of the droplet lateral to the directionof propagation. Fig. 5.5(b) shows the deformation of an identical two-dimensional incompressible dielectric droplet subjected to an s-polarizedGaussian beam with a peak electric field intensity of E0 = 0.98 GV/m.Each electrodynamic formalism predicts similar elongation of the dropletlongitudinal to the direction of propagation.(a) p, r = ℵ2, µr = 1 (b) s, r = ℵ2, µr = 1Figure 5.5: [Color online] Dielectric droplet deformation simulated usingPLBM. Given an s-polarized Hermite-Gaussian beam, each formalism pre-dicts similar interface deformation after t = 157.5 ns [t=150]. “REF” refersto the droplet’s initial shape.795.1. Two-dimensional simulationFig. 5.6(a) shows the deformation of a two-dimensional weakly compress-ible (Ma < 0.3) magnetic droplet (µr = 1.779), represented as a cylindricalcross section (or circular interface), subjected to an s-polarized Gaussianbeam with a peak electric field intensity of E0 = 0.98 GV/m. While theAB, MN, EL, HLM and AMP electrodynamic formalisms predict similarelongation longitudinal to the direction of propagation, the CHU formal-ism predicts elongation lateral to the direction of propagation. Fig. 5.6(b)shows the deformation of an identical two-dimensional incompressible mag-netic droplet subjected to a p-polarized Gaussian beam with a peak electricfield intensity of E0 = 0.98 GV/m. Most electrodynamic formalisms predictsimilar elongation of the droplet longitudinal to the direction of propagation,except for the AMP formalism.(a) s, r = 1, µr = ℵ2 (b) p, r = 1, µr = ℵ2Figure 5.6: [Color online] Magnetic droplet deformation simulated usingPLBM after t = 157.5 ns [t=150].“REF” refers to the droplet’s initial shape.5.1.4 ConclusionIn either the case of a flat interface or a circular interface, similar de-formation of the fluid is predicted by multiple electrodynamic formalismsdespite varying the lattice stencil range, the kinetic collision operator, thegeometry of the fluid, the geometry of the simulation domain, and theboundary conditions utilized in the PLBM. Therefore, the simulations inthis section represent the generic behaviour of dielectric or magnetic liq-uids with an optical index of water (ℵ = 1.333) when subjected to opticalradiation pressure from either s- or p-polarized Hermite-Gaussian beams.805.2. Three-dimensional simulation5.2 Three-dimensional simulationIn two dimensions the Chu and Amperian/Lorentz electrodynamic for-malisms predict a concave interface deflection from p-polarized optical ra-diation pressure, and a convex interface deflection from s-polarized opticalradiation pressure. It is therefore necessary to simulate Chu’s formalism inthree dimensions to properly analyze its polarization-dependent predictionsof interface deformation. Three-dimensional simulation is also necessaryfor representing elliptically polarized Gaussian beams, for which linearly-polarized and circularly-polarized Gaussian beams are but a special case.Linearly-polarized Gaussian beams in three dimensions also exhibit cross-polarization, or the existence of electric field components orthogonal tothe nominal polarization of the Gaussian beam [79]. Therefore, a three-dimensional simulation can represent both s- and p-polarized electric fieldsdepending on an analyzed plane of incidence. Although the cross-polarizedcomponents of an electric field have less intensity, they are still present,necessary to satisfy Maxwell’s equations, and can contribute to materialdeformation. Therefore, three-dimensional simulations of optical radiationpressure are recommended for their consistency with experiment.5.2.1 MethodologyELBM and PLBM are simpler to implement in three dimensions than theYL-FDTDM and front-tracking, even with adaptations for surface tensionand viscosity [92, 93]. The only difference between a two-dimensional andthree-dimensional case for the ELBM and PLBM are the collision operatorand the lattice vector set. Therefore, the deformation of spherical dielec-tric droplets from both linearly-polarized and circularly-polarized Gaussianbeams are analyzed in three dimensions.The spherical liquid-vapour interface, which is modeled using an SRTPLBM with a Shan-Chen forcing scheme, is used to predict the deformationof a “droplet” of refractive index ℵ = 1.333 and diameter D = 1.5λ0 sub-jected to optical radiation pressure from linearly-polarized and circularly-polarized Gaussian beams. Its spherical shape represents the most stablegeometry for a microscopic liquid-vapour interface. When initialized in thecenter of the simulation domain, it can also be acted on by optical forcesin all directions while avoiding kinetic interactions with the boundaries,which represent an external system. In each three-dimensional simulation,the liquid phase of the droplet is representative of water at 70 ◦C, and815.2. Three-dimensional simulationhas a dimensional mass density ρ′liq. = 977 kg/m3 and kinematic viscosityν′= 1.0034×10−6 m2/s. However, the surface tension is distinct from water:the dimensionless Shan-Chen surface tension was measured as γSC = 0.235with determination coefficient R2 = 0.9647. Fig. 5.7 shows an illustrationof the three-dimensional simulation configuration.m( 1060,0,1060) n−0 0w =zkxz3 6 2liq.Liquid Phase : 1.779, 1.333, 977 kg/m ' 1.0034 10 m /s' ,r  −= = = = 3 6 2vap.V s'apor phase : 1.00, 1.00, 1.00, kg/m ' 1.08 0340.8 0, 1 /2 mr r  −= = = = = m(1060,0, 060) n1m(1060,0, 1060) n−m( 1060,0, n1060)− −(0,0,0)01.5D =PBCPBCPBCPBCm(0, 1060,1060) n−0 0w =zkyzm(0,1060, 060) n1m(0,1060, 1060) n−m(0, 1060, n1060)− −(0,0,0)01.5D =PBCPBCPBCPBCFigure 5.7: [Color online] The configuration used in a three-dimensional sim-ulation of optical radiation pressure on a spherical interface. The dimensionsof the simulation domain are 4λ0× 4λ0× 4λ0, bordered by a black line rep-resenting periodic boundary conditions (PBCs). The gray region representsthe dielectric sphere, which has a diameter of D = 1.5λ0 and a permittivityof r = 1.779. The white region is free space. The green dashed line isthe source condition, which is a Hermite-Gaussian beam with a curvaturerepresented by the thin dashed black lines, and a beam waist of w0 = λ0(not to scale) centred at the position (0,0,0).825.2. Three-dimensional simulation5.2.2 The linearly-polarized Hermite-Gaussian beamThe linearly-polarized Hermite-Gaussian beam is polarized in thexˆ direction, and propagates in the zˆ direction, causing a deformation of thespherical droplet which is predicted by multiple elecrodynamic formalismsin Fig. 5.8. The electrodynamic formalisms can be divided into three classesdepending on the deformation of the droplet, which is dictated by materialparameters and their description within the momentum flux density tensor:1. The AB and MN electrodynamic formalisms predict identical defor-mation, since time-averaged time-variation of the time-harmonic mo-mentum density is zero. These two formalisms predict that the dropletdeforms from a spherical shape to a smooth ellipsoid with a major axisaligned to the direction of propagation.2. The EL and HLM electrodynamic formalisms predict identical defor-mation, since the medium is an isotropic dielectric. Similar to theAB/MN class, the EL/HLM class predicts that the droplet deformsfrom a sphere to an ellipsoid with a major axis aligned to the di-rection of propagation. However, the EL/HLM class predict time-averaged force density distributed throughout the dielectric, causingsmaller spatially periodic ripples in the interface which are not pre-dicted by the AB/MN class. This maybe an effect of electrostrictionat the micro-scale.3. The CHU and AMP electrodynamic formalisms predict identical ellip-soidal elongation of the droplet in the yˆ direction, which is trans-verse to the direction of propagation. The transverse fluid defor-mation was similarly predicted by Mansuripur, who suggested thata dielectric medium would form a “pancake” in his analysis of theLorentz/Amperian time-averaged force density distribution [16].Each class of electrodynamic formalism predicts that the droplet’s centerof mass accelerates in the direction of propagation, with a distinct dimen-sionless momentum transfer on the order of 102 (Fig. 5.11). Each electro-dynamic formalism also predicts zero dimensionless momentum transfer inlateral directions ( 10-14), validating the symmetry, isotropy, and momentumconservation of the simulation.835.2. Three-dimensional simulation(a) AB/MN (b) AB/MN (c) AB/MN(d) EL (e) EL (f) EL(g) HLM (h) HLM (i) HLM(j) CHU/AMP (k) CHU/AMP (l) CHU/AMPFigure 5.8: Three-dimensional deformation of an incompressible sphericaldroplet from an xˆ-polarized Gaussian beam with a peak electric field inten-sity of E0 = 0.98 GV/m (power P0 = 568 W). The dimensional kinematicviscosity is ν′= 1.0034× 10−6 m2/s, and the isosurface is plotted at dimen-sional time t′= 157.5 ns [t = 150] after the beam is turned on. 845.2. Three-dimensional simulation5.2.3 The circularly-polarized Hermite-Gaussian beamThe circularly-polarized Hermite-Gaussian beam propagates in thezˆ direction, with a rotating transverse electric field vector, causing a de-formation of the spherical droplet which is predicted by multiple electrody-namic formalisms in Fig. 5.9. The electrodynamic formalisms can be dividedinto three classes depending on the deformation of the droplet, which is dic-tated by material parameters and their description within the momentumflux density tensor:1. The AB and MN electrodynamic formalism predict that the dropletdeforms from a spherical shape to a smooth ellipsoid shape similar tothe smooth ellipsoid caused by a linearly-polarized Hermite-Gaussianbeam.2. The EL and HLM electrodynamic formalisms predict that the dropletdeforms from a spherical shape to a rippled ellipsoid shape similarto the rippled ellipsoid shape caused by a linearly-polarized Hermite-Gaussian beam.3. The CHU and AMP electrodynamic formalisms identically predict anellipsoidal elongation in the zˆ direction, which is the direction of prop-agation. Deformation in the direction of propagation is distinct fromthe CHU/AMP prediction of the droplet’s deformation from a lin-early polarized Hermite-Gaussian beam. The electic-field polariza-tion dependent deformation predicted by CHU and AMP formalismsresults in polarization dependent momentum transfer to the kineticsub-system over time. Given a circularly polarized Hermite-Gaussianbeam, the CHU/AMP formalisms also predict significantly less defor-mation along the direction of propagation, zˆ, compared with alterna-tive electrodynamic formalisms. Polarization-dependent fluid defor-mation has yet to be observed in experiment.Each class of electrodynamic formalism predicts that the droplet’s center ofmass accelerates in the direction of propagation, with a distinct dimension-less momentum transfer on the order of 102 (Fig. 5.11). Each electrody-namic formalism also predicts zero momentum transfer in lateral directions( 10-14), validating the symmetry, isotropy, and momentum conservation ofthe simulation.855.2. Three-dimensional simulation(a) AB/MN (b) AB/MN (c) AB/MN(d) EL (e) EL (f) EL(g) HLM (h) HLM (i) HLM(j) CHU/AMP (k) CHU/AMP (l) CHU/AMPFigure 5.9: Three-dimensional deformation of an incompressible sphericaldroplet from a circularly-polarized Gaussian beam with a peak electric fieldintensity of E0 = 0.98 GV/m (power P0 = 568 W). The dimensional kine-matic viscosity is ν′= 1.0034× 10−6 m2/s, and the isosurface is plotted atdimensional time t′= 157.5 ns [t = 150] after the beam is turned on. 865.2. Three-dimensional simulation5.2.4 Kinetic momentum transfer and electric fieldpolarizationDue to their unique predictions of droplet deformation, each class of for-malisms predicts a unique net dimensionless kinetic momentum along thelongitudinal zˆ direction (pk,z =∑xGk,zδx3), and zero transfer of kineticmomentum in lateral directions xˆ or yˆ, as illustrated in Figs. 5.10-5.11.The lateral kinetic momentum, pk,x =∑xGk,xδx3 and pk,y =∑xGk,yδx3,is theoretically zero. The numerical value of each lateral component illus-trated in Figs. 5.10-5.11 is negligible (10−14) relative to the dimensionlesskinetic momentum in the direction of propagation, pk,z. While most for-malisms predict the same momentum transfer independent of the electricfield polarization (Fig. 5.11), the Chu and Amperian/Lorentz formalismspredict a momentum transfer that is dependent on the electric field polar-ization, as illustrated in Fig. 5.10.(a) CHU/AMP, lin. (b) CHU/AMP, circ.Figure 5.10: The CHU/AMP formalisms momentum transfer for an (a)linearly-polarized Gaussian beam, and (b) a circularly-polarized Gaussianbeam. Each case predicts a distinct momentum transfer from the distinctshape of the droplet (pk,z,lin. = 220 while pk,z,circ. = 200 at t = 150). There-fore, the Chu and Amperian/Lorentz formalisms are the only electric fieldpolarization-dependent formalisms. There is no momentum transferred or-thogonal to the direction of propagation, due to symmetry. Error is on theorder of 10−14, and blue overlaps red.875.2. Three-dimensional simulation(a) AB/MN, lin. (b) AB/MN, circ.(c) EL, lin. (d) EL, circ.(e) HLM, lin. (f) HLM, circ.Figure 5.11: Dimensionless kinetic momentum from [left] a linearly-polarizedelectric field and [right] a circularly-polarized electric field. AB, MN, EL,MN formalisms predict polarization independent longitudinal momentum.Transverse momentum is on the order of 10−14, and blue overlaps red. 885.2. Three-dimensional simulation5.2.5 ConclusionFigs. 5.8 and 5.9 plot the deformation of a spherical droplet subjectedto linearly-polarized and circularly-polarized Gaussian beams, respectively.Time-averaged force density causes most formalisms to predict that thedroplet deforms along the direction of propagation, independent of electric-field polarization. The similar deformation predicted by the AB, MN, ELand HLM electrodynamic formalisms is illustrated by the paraboloid shapeof the droplet after t = 150 lattice time units. However, the Chu and Ampe-rian/Lorentz (CHU/AMP) formalisms can predict transverse polarization-dependent deformation instead. In each plot, the colorbar axis shows thatthe Mach number is below 0.3 (Ma < 0.3), such that the fluid remains in-compressible. However, weak compressibility can still occur at this Machnumber in lattice-Boltzmann simulations. The Helmholtz (HLM) and Ein-stein & Laub (EL) formalisms predict greater deformation of the dropletthan Abraham (AB) and Minkowski (MN) formalisms. Regardless, momen-tum is conserved in each numerical simulation with a dimensionless erroron the order of magnitude of 10-14, while the dimensionless net momentumtransferred from the electrodynamic to kinetic sub-system is on the order ofmagnitude of 102.89Chapter 6ConclusionSimulations attempting to isolate a single electrodynamic formalism ben-efit from alternative lattice-Boltzmann numerical methods which offer sim-pler implementation, potential increases in accuracy and precision, and de-creased computational time for representing unsteady incompressible flow intwo and three dimensions. Previous simulations of light-matter interactionsinvestigated in either two or three dimensions were limited to either a lowdensity ratio [7, 8], or linear approximations of a pressure at the interface[28]. In this thesis, Hauser & Verhey’s electrodynamic lattice Boltzmannmethod is integrated with both Shan & Chen single-relaxation time and Yu& Fan multi-relaxation-time pseudopotential lattice-Boltzmann methods tomodel the acceleration of generic two- and three-dimensional multi-phasefluids subjected to optical radiation pressure [39, 91, 93, 94]. When com-pared with the Yee-lattice finite-difference time-domain method, the electro-dynamic lattice-Boltzmann method evaluates a time-averaged kinetic forcedensity with greater accuracy in a decreased computational time. Single-component multi-phase fluids analyzed using the pseudopotential lattice-Boltzmann method achieve high density ratio and momentum conservationin two and three dimensions, while enabling the analysis of light-matterinteractions at the mesoscopic scale. The sub-system interaction is in agree-ment with previous studies and conclusions regarding optical radiation pres-sure in either two or three dimensions [7, 8, 16, 38]. Simulation results canbe scaled as necessary to represent experiments at characteristic lengthswhere the Navier-Stokes equation remains valid, and conditions for numeri-cal stability and low Mach number are maintained. Additional modificationshave also been suggested throughout this thesis for representing relativistichydrodynamics or higher Mach number and Reynolds number flows. Themethodology presented in this thesis shows that multiple electrodynamicformalisms can predict identical fluid-fluid interface deformation. However,this thesis finds that multiple electrodynamic formalisms predict distinctdeformation of linear dielectric media, deformations that are in qualitativeagreement with former studies and conclusions [7, 8, 16]. Therefore, thisthesis satisfies the following objectives:90Chapter 6. Conclusion1. Implement an alternative numerical approach for evaluating electro-dynamic force density in uniform isotropic media that improves onformer approaches.2. Implement a multi-physics simulation of optical radiation pressure act-ing on unsteady mesoscopic multi-phase fluids in two and three dimen-sions.3. Use the multi-physics simulation of optical radiation pressure to in-vestigate two and three-dimensional systems capable of representingsimilarities or differences between electrodynamic formalism predic-tions in the near-rest (quasi-stationary) frame of reference.Suitable numerical approaches in the literature were investigated, and thesatisfaction of the former objectives lead to the following conclusions:1. Lattice-Boltzmann methods, which consider particle population den-sities as opposed to phenomenological fields, are more accurate andsimpler to implement for solving electromagnetic force density in twoand three dimensions. Such methods have the same implementationregardless of the number of dimensions, and are devoid of lattice asym-metries or numerical dispersion in free space. At the same space andtime resolution, the electrodynamic lattice-Boltzmann method is 11 to35 times more accurate and 4 to 8 times faster than the Yee-latticefinite-difference time-domain method for evaluating a time-averagedkinetic force density.2. The same lattice-Boltzmann algorithm can be applied to represent ki-netic, electrodynamic and thermodynamic sub-systems by changingthe collision operation, lattice stencil, equilibrium condition and inter-acting forces. It’s adaptability, portability and parallel operations areconvenient for multi-physics simulation.3. Lattice-Boltzmann methods provide similar predictions to former method-ology, but their stability enables the simulation of greater density-ratiomesoscopic fluids with natural interface break-up, vaporization and co-alescence.4. Given a linearly-polarized Hermite-Gaussian beam transmitted througha quasi-stationary incompressible droplet, each electrodynamic for-malism (besides the Chu and Amperian/Lorentz formalisms) predictthat a spherical droplet suspended in its own vapor will form into a91Chapter 6. Conclusionparaboloid with its major axis oriented along the direction of propaga-tion. Given a circularly-polarized Gaussian beam transmitted througha quasi-stationary incompressible droplet, each of the formalisms (in-cluding the Chu and Amperian/Lorentz formalisms) predict the for-mation of this paraboloid.5. The Minkowski and Abraham formalisms predict time-averaged sur-face force densities that cause identical quasi-stationary isotropic mag-netodielectric fluid interface deformation from time-harmonic electro-magnetic waves. There is no difference in the time-averaged forces andfluid deformation occurring from the Abraham or Minkowski momen-tum density.6. The Helmholtz and Einstein & Laub electrodynamic formalisms pre-dict time-averaged bulk force densities that cause identical quasi-stationaryisotropic magnetodielectric fluid interface deformation. Their predic-tion is distinct from the Abraham and Minkowski electrodynamic for-malisms.7. The Chu and Amperian/Lorentz formalisms predict identical forcedensity in a dielectric medium, but distinct force density in a mag-netic medium.8. The Chu and Amperian/Lorentz electrodynamic formalisms predictthat the deformation of a quasi-stationary isotropic magnetodielectricfluid interface depends on the polarization of the incident electric field,while alternative electrodynamic formalisms do not. Consequently, theChu and Amperian/Lorentz electrodynamic formalisms also predict apolarization-dependent momentum transfer in a dielectric fluid overtime.9. Electrodynamic formalisms can be classified according to their interac-tion with a deformable medium. The cause of the distinct deformationspredicted by each class of formalism may be attributed to the distinctformulation of Maxwell’s equations from which they are derived. Dif-ferences in the description of fields internal to a medium will requirefurther investigation.Besides Chu’s formalism, most other electrodynamic formalisms provide de-generate predictions of liquid-vapor interface deformation that are largelyindependent of the electric field polarization. Differences between the for-malisms in each class may only exist in cases of compressible fluids and926.1. Recommendationswith their photothermal expansion [23, 28, 53]. Developing a mesoscopicthermodynamic sub-system was beyond the scope of this thesis, which onlyconsiders the case of an isothermal fluid. Therefore, thermodynamic or pho-toacoustic effects were not investigated, but are recommended for futurestudy.The systems described in this thesis are also non-dimensional, such thatthey do not represent particular fluid-fluid optical radiation pressure experi-ments, and they may use fluid properties and dimensions that are difficult torealize in experiment. However, the method’s dimensionless nature may bescaled to a fluid system of arbitrary size when provided with the necessarycomputational resources.6.1 RecommendationsIn future studies, the incorporation of the thermodynamic sub-systemand dimensional parameters would be useful for investigating experimentsor applications of optical radiation pressure on general lossy multiphase andmulti-component fluid systems. If performed as an academic research study,the scope of the development of this method for a particular experimentis recommended at the doctorate or post-doctorate level, preferably for amulti-disciplinary research group with independent members developing theindependent systems based on their disciplines. Significant computationalpower may also be required to represent the electrodynamic sub-system foroptically large simulation domains, but the resources necessary will dependon the system. The necessary computational resources require further in-vestigation.93Bibliography[1] M. Antognozzi, C. R. Bermingham, R. L. Harniman, S. Simpson, J. Se-nior, R. Hayward, H. Hoerber, M. R. Dennis, A. Y. Bekshaev, K. Y.Bliokh, and F. Nori. Direct measurements of the extraordinary op-tical momentum and transverse spin-dependent force using a nano-cantilever. Nature Physics, 12(8):731–735, 2016. ISSN 17452481. doi:10.1038/nphys3732. → pages 1[2] David Kipping. Relativistic light sails. arXiv, 153(6):277, 2017. ISSN23318422. doi: 10.3847/1538-3881/aa729d. → pages 1, 6[3] Duncan Farquharson Gregory, Robert Leslie Ellis, William Thom-son Baron Kelvin, and Norman Macleod Ferrers. The Cambridgeand Dublin Mathematical Journal, volume 4. Macmillan, Barclay andMacmillan, 1849. → pages 1[4] James Clerk Maxwell. VIII. A dynamical theory of the electromag-netic field. Philosophical Transactions of the Royal Society of London,155:459–512, January 1865. ISSN 0370-2316 (print), 2053-9207 (elec-tronic). doi: https://doi.org/10.1098/rstl.1865.0008. → pages 1, 12[5] A. Shevchenko and B. J. Hoenders. Microscopic derivation of elec-tromagnetic force density in magnetic dielectric media. New Journalof Physics, 12, 2010. ISSN 13672630. doi: 10.1088/1367-2630/12/5/053020. → pages 1, 2, 3, 14, 15, 74[6] Hendrik Antoon Lorentz. Versuch einer theorie der electrischen undoptischen erscheinungen in bewegten ko¨rpern, ej brill, leiden 1895.The English translation: Attempt of a theory of electrical and opticalphenomena in moving bodies can be found in: http://en. wikisource.org/wiki/Translation: Attempt of a Theory of Electrical and OpticalPhenomena in Moving Bodies, 1895. → pages 1[7] Max Bethune-Waddell and Kenneth J Chau. Simulations of radiationpressure experiments narrow down the energy and momentum of light94Bibliographyin matter. Reports on Progress in Physics, 78(12):122401, dec 2015.ISSN 0034-4885. doi: 10.1088/0034-4885/78/12/122401. → pages 1,2, 3, 4, 5, 9, 28, 30, 31, 33, 36, 38, 47, 54, 55, 56, 57, 73, 74, 77, 90[8] Maximilien Bethune-waddell. Simulations of Radiation Pressure Ex-periments. Text, University of British Columbia, 2016. → pages 1, 2,3, 5, 8, 9, 14, 15, 31, 33, 36, 38, 54, 55, 73, 77, 90[9] Peter Lebedew. Untersuchungen u¨ber die Druckkra¨fte des Lichtes.Annalen der Physik, 311(11):433–458, 1901. → pages 1[10] J H Poynting. XXXIX. Radiation pressure. The London, Edinburgh,and Dublin Philosophical Magazine and Journal of Science, 9(52):393–406, 1905. → pages 1, 10, 14[11] A. Einstein. Zur Elektrodynamik bewegter Ko¨rper. Annalen derPhysik, 322(10):891–921, 1905. ISSN 00033804. doi: 10.1002/andp.19053221004. → pages 2[12] Hendrik Antoon Lorentz. Electromagnetic phenomena in a systemmoving with any velocity smaller than that of light. In Collected Pa-pers, pages 172–197. Springer, 1937. → pages 2[13] H Minkowski. Die grundgleichungen fu¨r die elektromagnetischenvorga¨nge in bewegten. Nachr. Ko¨nigl. Ges. Wiss. Go¨ettingen, 20:53111, 1908. → pages 2, 14, 15[14] Max Abraham. Zur Elektrodynamik bewegter Ko¨rper. Rendicontidel Circolo Matematico di Palermo, 28(1):1–28, dec 1909. ISSN 0009-725X. doi: 10.1007/BF03018208. → pages 2, 14, 15[15] A Einstein and J Laub. U¨ber die im elektromagnetischen Felde aufruhende Ko¨rper ausgeu¨bten ponderomotorischen Kra¨fte. Annalen derPhysik, 331(8):541–550, 1908. doi: 10.1002/andp.19083310807. →pages 2, 12, 14, 15[16] Masud Mansuripur, Armis R Zakharian, and Ewan M Wright.Electromagnetic-force distribution inside matter. Physical Review A,88(2):023826, aug 2013. ISSN 1050-2947. doi: 10.1103/PhysRevA.88.023826. → pages 2, 3, 12, 38, 56, 74, 83, 90[17] I. Brevik. Experiments in phenomenological electrodynamics and theelectromagnetic energy-momentum tensor. Physics Reports, 52(3):95Bibliography133–201, may 1979. ISSN 03701573. doi: 10.1016/0370-1573(79)90074-7. → pages 2[18] Iver Brevik, Prachi Parashar, and K. V. Shajesh. Casimir force formagnetodielectric media. Physical Review A, 98(3):032509, sep 2018.ISSN 2469-9926. doi: 10.1103/PhysRevA.98.032509. → pages 2[19] Robert M Fano, Richard B Adler, and Lan Jen Chu. Electromagneticfields, energy, and forces. Taylor & Francis, 1963. → pages 2, 3, 14,15[20] L. J. Chu, H. A. Haus, and P. Penfield. The Force Density in Po-larizable and Magnetizable Fluids. Proceedings of the IEEE, 54(7):920–935, 1966. ISSN 15582256. doi: 10.1109/PROC.1966.4938. →pages 2, 8, 9, 13, 14, 15, 18, 74[21] W. Shockley and R. P. James. ”try simplest cases” discovery of ”hid-den momentum” forces on ”magnetic currents”. Phys. Rev. Lett., 18:876–879, May 1967. doi: 10.1103/PhysRevLett.18.876. URL https://link.aps.org/doi/10.1103/PhysRevLett.18.876. → pages 3[22] A. Ashkin and J. M. Dziedzic. Radiation Pressure on a Free LiquidSurface. Physical Review Letters, 30(4):139–142, jan 1973. ISSN 0031-9007. doi: 10.1103/PhysRevLett.30.139. → pages 3, 6, 28[23] Tomazˇ Pozˇar, Jernej Lalosˇ, Alesˇ Babnik, Rok Petkovsˇek, MaxBethune-Waddell, Kenneth J. Chau, Gustavo V.B. Lukasievicz, andNelson G.C. Astrath. Isolated detection of elastic waves driven by themomentum of light. Nature Communications, 9(1):1–11, 2018. ISSN20411723. doi: 10.1038/s41467-018-05706-3. → pages 3, 4, 5, 9, 14,15, 31, 38, 54, 56, 73, 93[24] Alexis Casner and Jean-Pierre Delville. Giant Deformations of aLiquid-Liquid Interface Induced by the Optical Radiation Pressure.Physical Review Letters, 87(5):054503, jul 2001. ISSN 0031-9007. doi:10.1103/PhysRevLett.87.054503. → pages 3[25] Alexis Casner and Jean Pierre Delville. Laser-Induced HydrodynamicInstability of Fluid Interfaces. Physical Review Letters, 90(14):4, 2003.ISSN 10797114. doi: 10.1103/PhysRevLett.90.144503. → pages 5[26] Iver Brevik and Simen A˜ Ellingsen. Detection of the Abraham forcewith a succession of short optical pulses. Physical Review A - Atomic,96BibliographyMolecular, and Optical Physics, 86(2), 2012. ISSN 10941622. doi:10.1103/PhysRevA.86.025801. → pages 54[27] Iver Brevik. What do experiments in optics tell us about photonmomentum in media? arXiv preprint arXiv:1310.3684, 2013. → pages3, 54, 74, 78[28] Nelson G. C. Astrath, Luis C. Malacarne, Mauro L. Baesso, GustavoV. B. Lukasievicz, and Stephen E. Bialkowski. Unravelling the effectsof radiation forces in water. Nature Communications, 5(1):4363, sep2014. ISSN 2041-1723. doi: 10.1038/ncomms5363. → pages 3, 5, 6, 9,15, 54, 57, 73, 78, 90, 93[29] Li Zhang, Weilong She, Nan Peng, and Ulf Leonhardt. Experimentalevidence for Abraham pressure of light. New Journal of Physics, 17(5):053035, may 2015. ISSN 1367-2630. doi: 10.1088/1367-2630/17/5/053035. → pages 4, 55[30] Mikko Partanen, Teppo Ha¨yrynen, Jani Oksanen, and Jukka Tulkki.Photon mass drag and the momentum of light in a medium. PhysicalReview A, 95(6):063850, jun 2017. ISSN 2469-9926. doi: 10.1103/PhysRevA.95.063850. → pages 4, 8, 15, 18[31] Brandon A. Kemp and Cheyenne J. Sheppard. Modeling opticalmanipulation using the field-kinetic and canonical formulations ofelectrodynamics. In Kishan Dholakia and Gabriel C. Spalding, ed-itors, Optical Trapping and Optical Micromanipulation XVI, vol-ume 11083, page 84. SPIE, sep 2019. ISBN 9781510628595. doi:10.1117/12.2528209. → pages 4, 5, 8[32] N. L. Balazs. The Energy-Momentum Tensor of the ElectromagneticField inside Matter. Physical Review, 91(2):408–411, jul 1953. ISSN0031-899X. doi: 10.1103/PhysRev.91.408. → pages 4[33] Kenneth J Chau and Henri J Lezec. Revisiting the Balazs thoughtexperiment in the case of a left-handed material: electromagnetic-pulse-induced displacement of a dispersive, dissipative negative-indexslab. Optics Express, 20(9):10138, apr 2012. ISSN 1094-4087. doi:10.1364/OE.20.010138. → pages 4, 5, 31[34] Robert N.C. Pfeifer, Timo A. Nieminen, Norman R. Heckenberg, andHalina Rubinsztein-Dunlop. Colloquium: Momentum of an electro-magnetic wave in dielectric media. Reviews of Modern Physics, 79(4):97Bibliography1197–1216, 2007. ISSN 00346861. doi: 10.1103/RevModPhys.79.1197.→ pages 5[35] N. G.C. Astrath, L. C. Malacarne, P. R.B. Pedreira, A. C. Bento,M. L. Baesso, and J. Shen. Time-resolved thermal mirror for nanoscalesurface displacement detection in low absorbing solids. Applied PhysicsLetters, 91(19), 2007. ISSN 00036951. doi: 10.1063/1.2812542. →pages 5[36] Zhi Tao Li, Gao Jin Li, Hai Bo Huang, and Xi Yun Lu. Latticeboltzmann study of electrohydrodynamic drop deformation with largedensity ratio. International Journal of Modern Physics C, 22(7):729–744, 2011. ISSN 01291831. doi: 10.1142/S0129183111016580. → pages5, 6, 55, 56, 73[37] S. M. Hanasoge, S. Succi, and S. A. Orszag. Lattice Boltzmannmethod for electromagnetic wave propagation. Epl, 96(1), 2011. ISSN02955075. doi: 10.1209/0295-5075/96/14002. → pages 24, 27, 40, 44[38] Ji Young Moon, Se Bin Choi, Jung Shin Lee, Roger I. Tanner, andJoon Sang Lee. Numerical simulation of optical control for a softparticle in a microchannel. Physical Review E, 99(2):1–10, 2019. ISSN24700053. doi: 10.1103/PhysRevE.99.022607. → pages 6, 55, 56, 73,90[39] A. Hauser and J. L. Verhey. Comparison of the lattice-Boltzmannmodel with the finite-difference time-domain method for electrody-namics. Physical Review E, 99(3):033301, mar 2019. ISSN 2470-0045.doi: 10.1103/PhysRevE.99.033301. → pages 22, 24, 27, 38, 39, 40, 41,44, 45, 47, 48, 51, 56, 90, 112[40] Li Chen, Qinjun Kang, Yutong Mu, Ya Ling He, and Wen Quan Tao.A critical review of the pseudopotential multiphase lattice Boltzmannmodel: Methods and applications. International Journal of Heat andMass Transfer, 76:210–236, 2014. ISSN 00179310. doi: 10.1016/j.ijheatmasstransfer.2014.04.032. → pages 5, 23, 55, 56, 58, 59, 61[41] Sauro Succi. The lattice Boltzmann equation: for fluid dynamics andbeyond. Oxford university press, 2001. → pages 5, 61[42] R Robert Nourgaliev, Truc-Nam Dinh, Theo G Theofanous, andD Joseph. The lattice boltzmann equation method: theoretical inter-98Bibliographypretation, numerics and implications. International Journal of Multi-phase Flow, 29(1):117–169, 2003. → pages 61[43] Bilen Emek Abali. An accurate finite element method for the numer-ical solution of isothermal and incompressible flow of viscous fluid.Fluids, 4(1), 2019. ISSN 23115521. doi: 10.3390/fluids4010005. →pages 5, 55[44] Sebastian Geller, Manfred Krafczyk, Jonas To¨lke, Stefan Turek, andJaroslav Hron. Benchmark computations based on lattice-boltzmann,finite element and finite volume methods for laminar flows. Computers& fluids, 35(8-9):888–897, 2006. → pages 5[45] Erlend Magnus Viggen. Acoustic equations of state for simple latticeboltzmann velocity sets. Physical Review E, 90(1):013310, 2014. →pages 5, 19, 20, 21, 22[46] I. V. Karlin, F. Bo¨sch, and S. S. Chikatamarla. Gibbs’ principle for thelattice-kinetic theory of fluid dynamics. Physical Review E - Statistical,Nonlinear, and Soft Matter Physics, 90(3), 2014. ISSN 15502376. doi:10.1103/PhysRevE.90.031302. → pages 17, 25, 55, 61, 65, 66[47] Huidan Yu and Kaihua Zhao. Lattice boltzmann method for com-pressible flows with high mach numbers. Physical Review E, 61(4):3867, 2000. → pages[48] Abolfazl Toghaniyan, Majid Zarringhalam, Omid Ali, Gholamreza Ah-madi Sheikh Shabani, and Davood Toghraie. Application of latticeBoltzmann method and spinodal decomposition phenomenon for sim-ulating two-phase thermal flows. Physica A, 509:673–689, 2018. doi:10.1016/j.physa.2018.06.030. → pages[49] M. Mendoza and J. D. Mun˜oz. Three-dimensional lattice Boltzmannmodel for electrodynamics. Physical Review E - Statistical, Nonlinear,and Soft Matter Physics, 82(5):1–13, 2010. ISSN 15393755. doi: 10.1103/PhysRevE.82.056708. → pages 18, 27, 40[50] Paul Romatschke. Relativistic (lattice) Boltzmann equation withnonideal equation of state. Physical Review D - Particles, Fields,Gravitation and Cosmology, 85(6):1–10, 2012. ISSN 15507998. doi:10.1103/PhysRevD.85.065012. → pages 5, 6, 1899Bibliography[51] A. N. Grigorenko, N. W. Roberts, M. R. Dickinson, and Y. Zhang.Nanometric optical tweezers based on nanostructured substrates. Na-ture Photonics, 2(6):365–370, 2008. ISSN 17494885. doi: 10.1038/nphoton.2008.78. → pages 6[52] A. Ashkin and J. M. Dziedzic. Internal cell manipulation using infraredlaser traps. Proceedings of the National Academy of Sciences of theUnited States of America, 86(20):7914–7918, 1989. ISSN 00278424.doi: 10.1073/pnas.86.20.7914. → pages 6, 57[53] A. Ashkin. Optical trapping and manipulation of small neutral parti-cles using lasers. Optics InfoBase Conference Papers, 94(May):4853–4860, 2010. ISSN 21622701. doi: 10.1364/fio.2010.stuc1. → pages 6,93[54] Dongmin Wang, Ping Cheng, and Xiaojun Quan. Photothermalnanobubble nucleation on a plasmonic nanoparticle: A 3D latticeBoltzmann simulation. International Journal of Heat and MassTransfer, 140:786–797, 2019. ISSN 00179310. doi: 10.1016/j.ijheatmasstransfer.2019.05.096. → pages 6, 66[55] Atsuo Murakami. Dynamics of Laser-Induced Bubble in PressurizedLiquid Nitrogen. JSME International Journal, 39(2):257–263, 1996.→ pages 6[56] Y. Tomita, M. Tsubota, K. Nagane, and N. An-naka. Behavior oflaser-induced cavitation bubbles in liquid nitrogen. Journal of AppliedPhysics, 88(10):5993–6001, nov 2000. ISSN 0021-8979. doi: 10.1063/1.1320028. → pages 6[57] Y. Tomita, M. Tsubota, and N. An-naka. Energy evaluation of cav-itation bubble generation and shock wave emission by laser focusingin liquid nitrogen. Journal of Applied Physics, 93(5):3039–3048, mar2003. ISSN 0021-8979. doi: 10.1063/1.1542669. → pages 6[58] M. Mendoza, B. M. Boghosian, H. J. Herrmann, and S. Succi. Fast lat-tice Boltzmann solver for relativistic hydrodynamics. Physical ReviewLetters, 105(1):2–5, 2010. ISSN 00319007. doi: 10.1103/PhysRevLett.105.014502. → pages 6, 55[59] B. U. Felderhof. Relativistic hydrodynamics of magnetic and di-electric fluids in interaction with the electromagnetic field. Journal100Bibliographyof Chemical Physics, 120(8):3598–3603, 2004. ISSN 00219606. doi:10.1063/1.1642612. → pages 6[60] Re´mi Soulard, Mark N. Quinn, Toshiki Tajima, and Ge´rard Mourou.ICAN: A novel laser architecture for space debris removal. Acta As-tronautica, 105(1):192–200, 2014. ISSN 00945765. doi: 10.1016/j.actaastro.2014.09.004. → pages 6[61] Vitaly V. Datsyuk and Oleg R. Pavlyniuk. Maxwell stress on a smalldielectric sphere in a dielectric. Physical Review A - Atomic, Molec-ular, and Optical Physics, 91(2):1–6, 2015. ISSN 10941622. doi:10.1103/PhysRevA.91.023826. → pages 10[62] Jean-Pierre Berenger. A perfectly matched layer for the absorptionof electromagnetic waves. Journal of Computational Physics, 114(2):185–200, 1994. ISSN 0021-9991. doi: 10.1006/jcph.1994.1159. →pages 12, 37[63] Masud Mansuripur. The force law of classical electrodynamics:Lorentz versus einstein and laub. In Optical Trapping and OpticalMicromanipulation X, volume 8810, page 88100K. International Soci-ety for Optics and Photonics, 2013. → pages 12[64] Cheyenne J. Sheppard and Brandon A. Kemp. Kinetic-energy-momentum tensor in electrodynamics. Physical Review A, 93(1):1–13,2016. ISSN 24699934. doi: 10.1103/PhysRevA.93.013855. → pages15[65] James Clerk Maxwell. Iv. on the dynamical theory of gases. Philosoph-ical transactions of the Royal Society of London, (157):49–88, 1867. →pages 16, 17, 21, 22, 65[66] Ludwig Boltzmann. Weitere studien u¨ber das wa¨rmegleichgewichtunter gasmoleku¨len. In Kinetische Theorie II, pages 115–225. Springer,1970. → pages 17, 18[67] Xiaowen Shan and Since Rothman. Lattice Boltzmann model for sim-ulating flows with multiple phases and components. Physical ReviewE - Statistical, Nonlinear, and Soft Matter Physics, 47(3):1815–1820,1993. → pages 17, 55, 56, 58[68] George Euge`ne Uhlenbeck, George W Ford, Gordon William Ford, andElliot Waters Montroll. Lectures in statistical mechanics, volume 1.Proquest/Csa Journal Division, 1963. → pages 17, 18, 19101Bibliography[69] Sydney Chapman, Thomas George Cowling, and D Burnett. TheMathematical Theory of Non-uniform Gases: An Account of KineticTheory of Viscosity, Thermal Conduction, and Diffusion. CambridgeUniversity Press, 1952. → pages 19[70] Guy R McNamara and Gianluigi Zanetti. Use of the boltzmann equa-tion to simulate lattice-gas automata. Physical review letters, 61(20):2332, 1988. → pages 20[71] FJ Higuera, S Succi, and R Benzi. Lattice gas dynamics with enhancedcollisions. EPL (Europhysics Letters), 9(4):345, 1989. → pages 24, 25,58[72] F J Higuera and Javier Jime´nez. Boltzmann approach to lattice gassimulations. EPL (Europhysics Letters), 9(7):663, 1989. → pages 20[73] Goncalo Silva and Viriato Semiao. Truncation errors and the rota-tional invariance of three-dimensional lattice models in the latticeboltzmann method. Journal of Computational Physics, 269:259–279,2014. → pages 20, 21, 22[74] Irina Ginzburg, Dominique d’Humie`res, and Alexander Kuzmin. Opti-mal stability of advection-diffusion lattice boltzmann models with tworelaxation times for positive/negative equilibrium. Journal of Statis-tical Physics, 139(6):1090–1143, 2010. → pages 22[75] Rashmi Ramadugu, Sumesh P Thampi, Ronojoy Adhikari, SauroSucci, and Santosh Ansumali. Lattice operators from discrete hy-drodynamics. arXiv preprint arXiv:1208.1009, 2012. → pages 23[76] Xiaoyi He and Li-Shi Luo. Theory of the lattice boltzmann method:From the boltzmann equation to the lattice boltzmann equation. Phys-ical review E, 56(6):6811, 1997. → pages 24, 26, 61, 65, 66[77] P. L. Bhatnagar, E. P. Gross, and M. Krook. A model for collisionprocesses in gases. I. Small amplitude processes in charged and neutralone-component systems. Physical Review, 94(3):511–525, 1954. ISSN0031899X. doi: 10.1103/PhysRev.94.511. → pages 25, 61[78] A. Hauser and J. L. Verhey. Stable lattice Boltzmann model forMaxwell equations in media. Physical Review E, 96(6):063306, dec2017. ISSN 2470-0045. doi: 10.1103/PhysRevE.96.063306. → pages27, 40, 41, 44102Bibliography[79] W. L. Erikson and Surendra Singh. Polarization properties of Maxwell-Gaussian laser beams. Physical Review E, 49(6):5778–5786, jun 1994.ISSN 1063-651X. doi: 10.1103/PhysRevE.49.5778. → pages 28, 29,30, 56, 81, 112[80] A. Taflove and S. C. Hagness. Computational Electrodynamics: TheFinite-Difference Time-Domain Method. Artech House, 3rd edition,2005. → pages 31, 37, 38, 54[81] Masud Mansuripur. Electromagnetic force and torque in ponderablemedia. Optics Express, 16(19):14821, 2008. ISSN 1094-4087. doi:10.1364/oe.16.014821. → pages 38[82] Yu Cheng, Guangzhi Chen, Xiang-Hua Wang, and Shunchuan Yang.Investigation of numerical dispersion with time step of the FDTDmethods: Avoiding erroneous conclusions. CoRR, abs/2001.10721,2020. URL https://arxiv.org/abs/2001.10721. → pages 38[83] Zhenhai Shao, Zhongxiang Shen, Qiuyang He, and Guowei Wei. Ageneralized higher order finite-difference time-domain method and itsapplication in guided-wave problems. IEEE transactions on microwavetheory and techniques, 51(3):856–861, 2003. → pages 38[84] Zhongqiang Xie, Chi-Hou Chan, and Bo Zhang. An explicit fourth-order staggered finite-difference time-domain method for maxwell’sequations. Journal of Computational and Applied Mathematics, 147(1):75–98, 2002. → pages 38[85] Costas D Sarris. Extending the stability limit of the FDTD methodwith spatial filtering. IEEE Microwave and Wireless Components Let-ters, 21(4):176–178, 2011. ISSN 15311309. doi: 10.1109/LMWC.2011.2105467. → pages 38[86] Xihao Li, Costas D. Sarris, and Piero Triverio. Overcoming thefdtd stability limit via model order reduction and eigenvalue pertur-bation. In 2014 IEEE MTT-S International Microwave Symposium(IMS2014), pages 1–3, 2014. doi: 10.1109/MWSYM.2014.6848408. →pages 38[87] Yaxun Liu and Costas D. Sarris. Dynamically adaptive mesh refine-ment FDTD: A stable and efficient technique for time-domain simu-lations. Antem/URSI 2006 - 12th International Symposium on An-103Bibliographytenna Technology and Applied Electromagnetics and Canadian RadioSciences Conference, Proceedings, pages 3–6, 2006. → pages 47[88] G. Tryggvason, B. Bunner, A. Esmaeeli, D. Juric, N. Al-Rawahi,W. Tauber, J. Han, S. Nas, and Y. J. Jan. A Front-TrackingMethod for the Computations of Multiphase Flow. Journal of Com-putational Physics, 169(2):708–759, 2001. ISSN 00219991. doi:10.1006/jcph.2001.6726. → pages 54[89] Nicolo Frapolli, Shyam S Chikatamarla, and Iliya V Karlin. Entropiclattice boltzmann model for compressible flows. Physical Review E, 92(6):061301, 2015. → pages 55[90] Michael R. Swift, E. Orlandini, W. R. Osborn, and J. M. Yeomans.Lattice boltzmann simulations of liquid-gas and binary fluid systems.Phys. Rev. E, 54:5041–5052, Nov 1996. doi: 10.1103/PhysRevE.54.5041. → pages 55[91] Peng Yuan and Laura Schaefer. Equations of state in a lattice Boltz-mann model. Physics of Fluids, 18(4), 2006. ISSN 10706631. doi:10.1063/1.2187070. → pages 55, 58, 59, 60, 61, 90[92] Daniel Lycett-Brown and Kai H Luo. Improved forcing scheme inpseudopotential lattice Boltzmann methods for multiphase flow at ar-bitrarily high density ratios. Physical Review E - Statistical, Non-linear, and Soft Matter Physics, 91(2), 2015. ISSN 15502376. doi:10.1103/PhysRevE.91.023305. → pages 55, 58, 59, 61, 62, 81[93] Zhao Yu and Liang Shih Fan. Multirelaxation-time interaction-potential-based lattice Boltzmann model for two-phase flow. Physi-cal Review E - Statistical, Nonlinear, and Soft Matter Physics, 82(4):1–14, 2010. ISSN 15393755. doi: 10.1103/PhysRevE.82.046708. →pages 59, 61, 62, 63, 64, 65, 81, 90[94] W Yu and R Mittra. Advanced FDTD Methods: Parallelization, Ac-celeration, and Engineering Applications, 2011. → pages 55, 90[95] James Clerk-Maxwell. On the dynamical evidence of the molecularconstitution of bodies. Nature, 11(279):357–359, 1875. → pages 56[96] GK Batchelor. Fluid mechanics. by ld landau and em lifshitz. 2nd en-glish edition. pergamon press, 1987. 539 pp.£ 45 or 29.50 (paperback).Journal of Fluid Mechanics, 205:593–594, 1989. → pages 57104Bibliography[97] Zhaoli Guo, Chuguang Zheng, and Baochang Shi. Discrete latticeeffects on the forcing term in the lattice boltzmann method. Phys.Rev. E, 65:046308, Apr 2002. doi: 10.1103/PhysRevE.65.046308. →pages 58, 61, 62[98] Cheng Peng, Luis F. Ayala, Zhicheng Wang, and Orlando M. Ayala.Attainment of rigorous thermodynamic consistency and surface ten-sion in single-component pseudopotential lattice Boltzmann modelsvia a customized equation of state. Physical Review E, 101(6):63309,2020. ISSN 2470-0045. doi: 10.1103/physreve.101.063309. → pages59[99] Feifei Qin, Ali Mazloomi Moqaddam, Qinjun Kang, DominiqueDerome, and Jan Carmeliet. Entropic multiple-relaxation-time mul-tirange pseudopotential lattice Boltzmann model for two-phase flow.Physics of Fluids, 032104(30):1–11, 2018. doi: 10.1063/1.5016965. →pages 59, 61, 65, 66[100] Zhaoli Guo, Chuguang Zheng, and T. S. Zhao. A lattice BGK schemewith general propagation. Journal of Scientific Computing, 16(4):569–585, 2001. ISSN 08857474. doi: 10.1023/A:1013280900427. → pages59[101] Kannan N. Premnath and John Abraham. Three-dimensional multi-relaxation time (MRT) lattice-Boltzmann models for multiphase flow.Journal of Computational Physics, 224(2):539–559, 2007. ISSN10902716. doi: 10.1016/j.jcp.2006.10.023. → pages 62, 63[102] Q. Li, K. H. Luo, and X. J. Li. Lattice boltzmann modeling of mul-tiphase flows at large density ratio with an improved pseudopotentialmodel. Phys. Rev. E, 87:053301, May 2013. doi: 10.1103/PhysRevE.87.053301. → pages 62, 64, 65[103] A. Xu, T.S. Zhao, L. An, and L. Shi. A three-dimensional pseudo-potential-based lattice boltzmann model for multiphase flows withlarge density ratio and variable surface tension. International Jour-nal of Heat and Fluid Flow, 56:261–271, 2015. ISSN 0142-727X. doi:https://doi.org/10.1016/j.ijheatfluidflow.2015.08.001. → pages 62, 65[104] Sami Ammar, Guillaume Pernaudat, and Jean-Yves Tre´panier. A mul-tiphase three-dimensional multi-relaxation time (mrt) lattice boltz-mann model with surface tension adjustment. Journal of Compu-105Bibliographytational Physics, 343:73–91, 2017. ISSN 0021-9991. doi: https://doi.org/10.1016/j.jcp.2017.04.045. → pages 62, 65[105] Pierre Lallemand and Li-Shi Luo. Theory of the lattice boltzmannmethod: Dispersion, dissipation, isotropy, galilean invariance, and sta-bility. Phys. Rev. E, 61:6546–6562, Jun 2000. doi: 10.1103/PhysRevE.61.6546. URL https://link.aps.org/doi/10.1103/PhysRevE.61.6546. → pages 62, 63[106] Pierre Lallemand and Li-Shi Luo. Theory of the lattice boltzmannmethod: Acoustic and thermal properties in two and three dimensions.Phys. Rev. E, 68:036706, Sep 2003. doi: 10.1103/PhysRevE.68.036706.→ pages 63[107] Qing Li and K. H. Luo. Achieving tunable surface tension in thepseudopotential lattice boltzmann modeling of multiphase flows. Phys.Rev. E, 88:053307, Nov 2013. doi: 10.1103/PhysRevE.88.053307. →pages 65[108] F. Bo¨sch, B. Dorschner, and I. Karlin. Entropic multi-relaxation free-energy lattice boltzmann model for two-phase flows. EPL (EurophysicsLetters), 122(1):14002, apr 2018. doi: 10.1209/0295-5075/122/14002.→ pages 65, 66[109] Kai H Luo, Linlin Fei, and Geng Wang. A unified lattice boltzmannmodel and application to multiphase flows. Philosophical Transactionsof the Royal Society A, 379(2208):20200397, 2021. → pages 65, 66[110] Robert S Maier, Robert S Bernard, and Daryl W Grunau. Boundaryconditions for the lattice boltzmann method. Physics of Fluids, 8(7):1788–1801, 1996. → pages 67, 68[111] Anthony JC Ladd. Numerical simulations of particulate suspensionsvia a discretized boltzmann equation. part 1. theoretical foundation.Journal of fluid mechanics, 271:285–309, 1994. → pages 68[112] Qisu Zou and Xiaoyi He. On pressure and velocity boundary conditionsfor the lattice boltzmann bgk model. Physics of fluids, 9(6):1591–1598,1997. → pages 69[113] Attila Ma´rkus and Ga´bor Ha´zi. Simulation of evaporation by an ex-tension of the pseudopotential lattice boltzmann method: A quanti-tative analysis. Phys. Rev. E, 83:046705, Apr 2011. doi: 10.1103/106BibliographyPhysRevE.83.046705. URL https://link.aps.org/doi/10.1103/PhysRevE.83.046705. → pages 69[114] N. G.C. Astrath, G. V.B. Lukasievicz, L. C. Malacarne, and S. E.Bialkowski. Surface deformation effects induced by radiation pressureand electrostriction forces in dielectric solids. Applied Physics Letters,102(23), 2013. ISSN 00036951. doi: 10.1063/1.4809924. → pages 73[115] P. Penfield and H. A. Haus. Electrodynamics of Moving Media.Proceedings of the IEEE, 53(4):422, 1965. ISSN 15582256. doi:10.1109/PROC.1965.3795. → pages 74[116] Brandon A. Kemp. Macroscopic Theory of Optical Momentum.Progress in Optics, 60:437–488, 2015. ISSN 00796638. doi: 10.1016/bs.po.2015.02.005. → pages 74[117] Yangang Liu and Peter H Daum. Relationship of refractive index tomass density and self-consistency of mixing rules for multicomponentmixtures like ambient aerosols. Journal of Aerosol Science, 39(11):974–986, 2008. → pages 78[118] I. Brevik. Fluids in Electric and Magnetic Fields - Pressure Variationand Stability. Canadian journal of physics, 60(4):449–455, 1981. ISSN00084204. doi: 10.1139/p82-063. → pages 78107Appendix AAnalytical methodsIn this thesis, a variety of mathematical operations and analytical ap-proaches are used, including matrix operations, Taylor series expansion, andChapman-Enskog analysis. Special functions such as Hermite-polynomialsand/or Lyapunov functions are required to derive general equations or verifythe stability of numerical models. This section briefly reviews the mathe-matical operations and analytical methods used in this thesis.A.1 Coordinate systemsThis thesis uses a Cartesian coordinate system with absolute space andtime which are invariant with respect to one another. The coordinate sys-tem is also at rest. While the position vector, x, depends on the origin’slocation, (x0, y0, z0), the displacement vector, r, does not. Fig. A.1 showsan illustration of position and displacement vectors.Position vectorx = {x1xˆ1 + x2xˆ2 + x3xˆ3}= {xxˆ + yyˆ + zzˆ} . (A.1)Displacement vectorr = (x− x0)xˆ + (y − y0)yˆ + (z − z0)zˆ. (A.2)108A.2. Matrix operationsx 'xr 'r0 00 '(a) (b)Figure A.1: There are two different coordinate systems (a) and (b), wherethe coordinate system in (b) has been translated to (x0, y0, z0) from thecoordinate system in (a). Given the coordinate translation, x 6= x′, butr = r′.A.2 Matrix operationsMatrix operations used in this thesis include the divergence, the dyadicproduct, and the scalar product. Special matrices include the Kronecker-Delta or unit matrix, and the Levi-Civita matrix.Matrices can be used to describe a set of scalar components or orthogonalvectors in a multi-dimensional space. They may or may not describe physicalquantities, being only a method for organizing and performing operationswith sequences of numbers.Tensors are multi-dimensional matrices describing a linear relationshipbetween physical variables related to a vector space. Tensors are requiredto describe a state of stress at a point in two or three dimensions. Tensorelements can be represented with dimensional index α for a row, and dimen-sional index β for a column. Written in matrix form, a three-dimensionaltensor isTαβ ∈ T, [T] =T11 T12 T13T21 T22 T23T31 T32 T33 . (A.3)109A.2. Matrix operationsDivergence represents the outward flux of a vector or a tensor at a givenpoint. In the case of a vector, the divergence is a scalar function,∇ ·T = ∂∂ x1T1 +∂∂ x2T2 +∂∂ x3T3. (A.4)The divergence of a tensor is a vector that has both magnitude and direction,∇ ·T = ∂∂ xβTαβxˆα = ∂∂ x1T11 + ∂∂ x2T12 + ∂∂ x3T13∂∂ x1T21 + ∂∂ x2T22 + ∂∂ x3T23∂∂ x1T31 +∂∂ x2T32 +∂∂ x3T33 , (A.5)where xˆα represents a unit vector, pointed in the direction of an increasingdimension xα. In this thesis, x1 = x, x2 = y, and x3 = z. The divergenceof a tensor is necessary to describe the kinetic force density from opticalradiation pressure.Dyadic products of two equal-size tensors represent their outer product,T2= TT = T : T = [T][T]ᵀ, (A.6)where ᵀ represents a transpose operation, that swaps columns and rows of atensor. The same operation can be applied to equal-size vectors or matrices.Scalar products of two equal-size tensors represent their inner product,T ·T = [T]ᵀ[T], (A.7)where ᵀ represents a transpose operation, that swaps the columns and rowsof a tensor. The same operation can be applied to equal-size vectors ormatrices.The Kronecker-Delta or unit matrix can be represented as a squarediagonal matrix with each diagonal element equal to one,δαβ ={1, α = β,0, α 6= β, [I] =1 0 00 1 00 0 1 . (A.8)The unity matrix may also be referred to as the identity matrix.110A.3. Polynomials and functionsThe Levi-Civita matrix maps elements representing the sign for the per-mutation of dimensional indices, such that they have total antisymmetry.Given two dimensional indices, α and β,εαβ =+1, {α, β} = {1, 2},−1, {α, β} = {2, 1},0, α = β,(A.9)while given three dimensional indices, α, β and γ,εαβγ =+1, {α, β, γ} = {1, 2, 3}, {2, 3, 1}, {3, 1, 2},−1, {α, β, γ} = {3, 2, 1}, {1, 3, 2}, {2, 1, 3},0, α = β, β = γ, γ = α.(A.10)The Levi-Civita matrices in two and three dimensions,[Υ2D] =[0 +1−1 0]and [Υ3D] = 0 +1 −1−1 0 +1+1 −1 0 , (A.11)resemble the co-factor matrices used in the expansion of a determinant.Therefore, the Levi-Civita symbol or matrix is useful for mapping the cross-product between multiple vectors represented within tensors, such as theelectrodynamic lattice-Boltzmann method.A.3 Polynomials and functionsApproximations of functions can be represented using Taylor series poly-nomials and Hermite polynomials as detailed in this section. Lyapunovfunctions are also defined.Taylor series polynomials represent the sum of a function’s derivativesat a point, which converges to the function itself for an infinite number ofterms,f(x) =∞∑m=0f (m)(a)m!(x− a)m. (A.12)Taylor series having a finite number of terms may be used to derivecentral difference equations. Given a finite number of elements, Nm, a Taylor111A.3. Polynomials and functionsseries is a polynomial representation of a real function having an order ofaccuracy proportional to the number of its terms, O((x− a)Nm+1),f(x) =Nm∑m=0f (m)(a)m!(x− a)m +O((x− a)Nm+1). (A.13)Hermite polynomials approximate Gaussian distribution functions. Thesequence of Hermite polynomials are described byH(n)(x) = (−1)n 1w(x)dndxnw(x), w(x) =1√2pie−x22 (A.14)The first six Hermite polynomials are listed in Table A.1.Table A.1: Hermite polynomialsEven OddH(0)(x) = 1 H(1)(x) = xH(2)(x) = x2 − 1 H(3)(x) = x3 − 3xH(4)(x) = x4 − 6x2 + 3 H(5)(x) = x5 − 10x3 + 15xHermite polynomials can also be extended to d spatial dimensions,H(n)(x) = (−1)n 1w(x)∇(n)w(x), w(x) = 1(2pi)d2e−x2 , (A.15)wherein each Hermite polynomial is orthogonal with respect to w(x). Multi-dimensional Hermite polynomials are used to derive the second-order dis-crete Maxwell-Boltzmann distribution function representing the equilibriumcondition of an ideal gas, in addition to the electric field amplitude profilefor a three-dimensional Gaussian beam [79].Lyapunov functions are scalar functions used to prove the stable equi-librium of an ordinary differential equation. The existence of Lyapunovfunctions can be a necessary and adequate condition for stability. The Lya-punov function for an autonomous dynamical system isg|Rn → Rn, dxdt= g(x), (A.16)with an equilibrium point at x = 0. Lyapunov functions were used toderive a stable electrodynamic lattice-Boltzmann method, and are sufficientto demonstrate its stability [39].112Appendix BMomentum flux densitytensors expandedIn this section, the common electrodynamic momentum flux density ten-sors are expanded into a single matrix form. This is required to derive theirdivergence, and program their central differences.B.1 Abraham’s momentum flux density tensorT =12−ExDx +DyEy +DzEz−HxBx +HyBy +HzBz ,−DxEy − BxHy − ExDy −HxBy, −DxEz − BxHz − ExDz −HxBz ,−DyEx − ByHx − EyDx −HyBx,ExDx −DyEy +DzEz+HxBx −HyBy +HzBz ,−DyEz − ByHz − EyDz −HyBz ,−DzEx − BzHx − EzDx −HzBx, −DzEy − BzHy − EzDy −HzBy,ExDx +DyEy −DzEz+HxBx +HyBy −HzBzB.2 Minkowski’s momentum flux density tensorT =− 12DxEx +12DyEy +12DzEz− 12BxHx +12ByHy +12BzHz ,−DxEy − BxHy, −DxEz − BxHz ,−DyEx − ByHx,12DxEx −12DyEy +12DzEz+12BxHx −12ByHy +12BzHz ,−DyEz − ByHz ,−DzEx − BzHx, −DzEy − BzHy,12DxEx +12DyEy −12DzEz+12BxHx +12ByHy −12BzHz113B.3. Einstein & Laub’s momentum flux density tensorB.3 Einstein & Laub’s momentum flux densitytensorT =02(Ex2+ Ey2+ Ez2)+µ02(Hx2+Hy2+Hz2)−DxEx − BxHx,−DxEy − BxHy, −DxEz − BxHz ,−DyEx − ByHx,02(Ex2+ Ey2+ Ez2)+µ02(Hx2+Hy2+Hz2)−DyEy − ByHy,−DyEz − ByHz ,−DzEx − BzHx, −DzEy − BzHy,02(Ex2+ Ey2+ Ez2)+µ02(Hx2+Hy2+Hz2)−DzEz − BzHzB.4 Chu’s momentum flux density tensorT =02(−Ex2 + Ey2 + Ez2)+µ02(−Hx2 +Hy2 +Hz2),−0ExEy − µ0HxHy, −0ExEz − µ0HxHz ,−0EyEx − µ0HyHx,02(Ex2 − Ey2 + Ez2)+µ02(Hx2 −Hy2 +Hz2),−0EyEz − µ0HyHz ,−0EzEx − µ0HzHx, −0EzEy − µ0HzHy,02(Ex2+ Ey2 − Ez2)+µ02(Hx2+Hy2 −Hz2)B.5 Ampere’s momentum flux density tensorT =02(−Ex2 + Ey2 + Ez2)+12µ0(−Bx2 + By2 + Bz2),−0ExEy − 1µ0 BxBy, −0ExEz −1µ0BxBz ,−0EyEx − 1µ0 ByBx,02(Ex2 − Ey2 + Ez2)+12µ0(Bx2 − By2 + Bz2),−0EyEz − 1µ0 ByBz ,−0EzEx − 1µ0 BzBx, −0EzEy −1µ0BzBy,02(Ex2+ Ey2 − Ez2)+12µ0(Bx2+ By2 − Bz2)114B.6. Helmholtz momentum flux density tensorB.6 Helmholtz momentum flux density tensorT =CM (Ex2+ Ey2+ Ez2)− E2x+µCM (Hx2+Hy2+Hz2)− µH2x,−ExEy − µHxHy, −ExEz − µHxHz ,−EyEx − µHyHx,CM (Ex2+ Ey2+ Ez2)− E2y+µCM (Hx2+Hy2+Hz2)− µH2y ,−EyEz − µHyHz ,−EzEx − µHzHx, −EzEy − µHzHy,CM (Ex2+ Ey2+ Ez2)− E2z+µCM (Hx2+Hy2+Hz2)− µH2zwhere CM and µCM are the permittivity and the permeability for a Clausius-Mossotti medium, while  and µ are the standard permittivity and perme-ability. Each are given byCM = 02 + 2r − r26,  = 0r, (B.1)µCM = µ02 + 2µr − µr26, µ = µ0µr. (B.2)B.7 Divergence of a stress-tensorThe divergence of a stress-tensor is used to calculate the force density ona medium. The row and column indices are represented using α = {1, 2, 3}and β = {1, 2, 3}. Substituting α = {x, y, z} and β = {x, y, z} is a commonpractise for representing the row and column indices, such that the stress-tensor is represented as,Tαβ ∈ T, [T] =T11 T12 T13T21 T22 T23T31 T32 T33 =Txx Txy TxzTyx Tyy TyzTzx Tzy Tzz . (B.3)Divergence of this stress-tensor can then be represented as∇ ·T = ∂∂ xβTαβeα =∂∂ xTxx +∂∂ yTxy +∂∂ xTxz∂∂ xTyx +∂∂ yTyy +∂∂ yTyz∂∂ xTzx +∂∂ yTzy +∂∂ zTzz . (B.4)The divergence operator applies to all dimensions in which there is spatialvariation of the electric and magnetic fields. In the case of a two-dimensionalxy-plane simulation, spatial variation only occurs in the xy-plane, such thatthe terms ∂∂ xTxz = 0,∂∂ yTyz = 0,∂∂ xTzx = 0,∂∂ yTzy = 0 and∂∂ zTzz = 0.Therefore, the divergence of the stress-tensor may be represented as115B.7. Divergence of a stress-tensor∇ ·T = ∂∂ xTxx + ∂∂ yTxy + 0∂∂ xTyx + ∂∂ yTyy + 00 + 0 + 0 , (B.5)which demonstrates zero variation along the z-coordinate. The three-dimensionalvector with no information contained in one of its dimensions may be rep-resented as a two-dimensional vector,(∇ ·T)xy =[∂∂ xTxx +∂∂ yTxy∂∂ xTyx +∂∂ yTyy]=[∆Tx∆Ty]. (B.6)Two examples of the divergence are provided below for the Minkowskiand Chu formalisms. These two formalisms are chosen for their opposinginterpretation of light-matter interactions.B.7.1 Two-dimensional Minkowski stress-tensor divergence∇ · T = 12 ∂∂x(−DxEx +DyEy +BzHz) + ∂∂y (−2DxEy)∂∂x(−2DyEx) + ∂∂y (DxEx −DyEy +BzHz) . (B.7)B.7.2 Two-dimensional Chu stress-tensor divergence∇ · T = 12 ∂∂x( ε02 (−Ex2 + Ey2) + µ02 Hz2) + ∂∂y (−ε0ExEy)∂∂x(−ε0EyEx) + ∂∂y ( ε02 (Ex2 − Ey2) + µ02 Hz2) (B.8)116Appendix CDimensional ConversionsTab. C.1 provides a full list of the dimensional conversions used in boththe electrodynamic and kinetic sub-systems.Table C.1: Dimensional conversion of system variablesVariable Name Dimensionless (x, t) Dimensional (x′, t′)Time (e) te t′e = te∆t/2Time (k) tk t′k = tk∆tPosition vector x x′= x∆xVelocity u u′= u(∆x/∆t)Mass density ρ ρ′= ρρ′cElectric field intensity E E′= EMagnetic field intensity H H′= H/η0Energy density We W′e = We0Poynting’s vector Se S′e = Se/η0Momentum density (e) Ge G′e = Ge/(c2η0)Momentum density (k) Gk G′k = Gk(ρ′c∆x/∆t)Momentum-flux tensor Te T′e = Te0Force density fk f′k = fk[ρ′c∆x/(∆t)2]117Appendix DProgram ImplementationsThis section contains several program implementations of various MAT-LAB functions used in the ELBM and PLBM representing electrodynamicand kinetic sub-systems. Note that the functions are provided for the caseof a two-dimensional simulation.D.1 Electrodynamic lattice-Boltzmann methodforce densityThe following script is an excerpt of the electrodynamic lattice-Boltzmannmethod’s force density using simple central difference formulas. It is moreconvenient to write in an expanded form.% Calculate a kinetic force density from an electrodynamic formalism% model: indicates the electrodynamic formalism with an acronymn% er: relative permittivity% ur: relative permeability% dtp: dimensional time unit% dxp: dimensional space unit% Gp: Momentum density at previous time step% NE: number of electric field componentsfunction [F,G] = force(model,nx,ny,er,ur,dtp,E,H,Gp,NE,dxp)e0=8.8541878128e-12; %free space permittivityu0=1.25663706e-6; %free space permeabilityc=299792458; %speed of lighteta0=sqrt(u0/e0); %free space impedanceH=H/(2*eta0); %dimensional conversionE=E/2; %dimensional conversion% Initialize arrays according to the electric field polarizationif NE>1F=zeros(size(E)); G=zeros(size(E));elseF=zeros(size(H)); G=zeros(size(H));endii=4:nx-4;jj=4:ny-4;% Momentum densityif NE>1G(jj,ii,1)=E(jj,ii,1).*H(jj,ii,1)/(c^2);118D.1. Electrodynamic lattice-Boltzmann method force densityG(jj,ii,2)=-E(jj,ii,2).*H(jj,ii,1)/(c^2);elseG(jj,ii,1)=-E(jj,ii,1).*H(jj,ii,1)/(c^2);G(jj,ii,2)=E(jj,ii,1).*H(jj,ii,2)/(c^2);endF=zeros(size(G));% Abraham force density, polarization dependentif strcmp(model, ’AB’)==1% for a TM modeif NE>1t11 = (1/2)*(e0*er.*(sum(E.^2,3)-2*E(:,:,2).^2)+u0*ur.*(sum(H.^2,3)));t22 = (1/2)*(e0*er.*(sum(E.^2,3)-2*E(:,:,1).^2)+u0*ur.*(sum(H.^2,3)));t12 = e0*(-er.*E(:,:,1).*E(:,:,2)); t21=t12;elset11 = (1/2)*(e0*er.*(sum(E.^2,3))+u0*ur.*(sum(H.^2,3)-2*H(:,:,2).^2));t22 = (1/2)*(e0*er.*(sum(E.^2,3))+u0*ur.*(sum(H.^2,3)-2*H(:,:,1).^2));t12 = u0*(-ur.*H(:,:,1).*H(:,:,2)); t21=t12;end% Einstein & Laub force density, polarization dependentelseif strcmp(model,’EL’)==1if NE>1t11 = (1/2)*(e0*sum(E.^2,3)+u0*sum(H.^2,3)) - e0*er.*E(:,:,2).^2;t12 = (-e0*er.*E(:,:,1).*E(:,:,2)); t21=t12;t22 = (1/2)*(e0*sum(E.^2,3)+u0*sum(H.^2,3)) - e0*er.*E(:,:,1).^2;elset11 = (1/2)*(e0*sum(E.^2,3)+u0*sum(H.^2,3)) - u0*ur.*H(:,:,2).^2;t12 = (-u0*ur.*H(:,:,1).*H(:,:,2)); t21=t12;t22 = (1/2)*(e0*sum(E.^2,3)+u0*sum(H.^2,3)) - u0*ur.*H(:,:,1).^2;end% Chu force density, polarization dependentelseif strcmp(model,’CHU’)==1% for a TM modeif NE>1t11 = e0*(1/2)*(sum(E.^2,3)-2*E(:,:,2).^2)+u0*(1/2)*(sum(H.^2,3));t12 = e0*(-E(:,:,1).*E(:,:,2)); t21=t12;t22 = e0*(1/2)*(sum(E.^2,3)-2*E(:,:,1).^2)+u0*(1/2)*(sum(H.^2,3));elset11 = e0*(1/2)*(sum(E.^2,3))+u0*(1/2)*(sum(H.^2,3)-2*H(:,:,2).^2);t12 = -u0*H(:,:,1).*H(:,:,2); t21=t12;t22 = e0*(1/2)*(sum(E.^2,3))+u0*(1/2)*(sum(H.^2,3)-2*H(:,:,1).^2);end%Helmholtz force density, polarization dependentelseif strcmp(model, ’HLM’)==1ecm = (2+2*er-er.^2)/6;ucm = (2+2*ur-ur.^2)/6;if NE>1t11 = e0*ecm.*sum(E.^2,3)-e0*er.*E(:,:,2).^2+u0*ucm.*H(:,:,1).^2;t12 = e0*(-er.*E(:,:,1).*E(:,:,2)); t21=t12;t22 = e0*ecm.*sum(E.^2,3)-e0*er.*E(:,:,1).^2+u0*ucm.*H(:,:,1).^2;elset11 = u0*ucm.*sum(H.^2,3)-u0*ur.*H(:,:,2).^2+e0*ecm.*(E(:,:,1).^2);t12 = u0*(-ur.*H(:,:,1).*H(:,:,2)); t21=t12;t22 = u0*ucm.*sum(H.^2,3)-u0*ur.*H(:,:,1).^2+e0*ecm.*(E(:,:,1).^2);119D.2. Shan-Chen force density with applicable equation of stateendendif strcmp(model, ’HLM’)==1 %Clausius-Mossoti medium parameterseucm = (2+2*er(jj,ii).*ur(jj,ii)-2*er(jj,ii)+ur(jj,ii))/3;for mm=1:length(ii)for nn=1:length(jj)G(nn,mm,1)=G(nn,mm,1)*eucm(nn,mm);G(nn,mm,2)=G(nn,mm,2)*eucm(nn,mm);endendenddTx=zeros(size(G)); dTy=zeros(size(G));dTx(jj,ii) = (t11(jj,ii+1)-t11(jj,ii-1))/2 ...+(t12(jj+1,ii)-t12(jj-1,ii))/2;dTy(jj,ii) = (t21(jj,ii+1)-t21(jj,ii-1))/2 ...+(t22(jj+1,ii)-t22(jj-1,ii))/2;for mm=jjfor nn=iiF(mm,nn,1) = - (dTx(mm,nn)/(dxp) + (G(mm,nn,1)-Gp(mm,nn,1))/(2*dtp));F(mm,nn,2) = - (dTy(mm,nn)/(dxp) + (G(mm,nn,2)-Gp(mm,nn,2))/(2*dtp));endendendD.2 Shan-Chen force density with applicableequation of stateModulation operators are used to distribute the velocity vectors across periodic boundaries.This enables the definition of the gradient operator in the Shan-Chen force density throughoutall lattice nodes in the simulation domain. However, the implementation for a multi-range latticewould have greater implementation complexity.% Shan-Chen force densityfunction [fx,fy]=force(nx,ny,cx,cy,rho,w,G,fx,fy)% Equation of state parameters: rhov/rhol = 1180a=1; b=4; R=1;Tc=0.3773*a/(b*R); T=0.531*Tc;for j=1:nyfor i=1:nxfxs=0.0;fys=0.0;for k=1:8newx=1+mod(i-1+cx(k)+nx,nx);newy=1+mod(j-1+cy(k)+ny,ny);r=rho(newx,newy);p = (r*R*T*(1+(b*r/4)+(b*r/4)^2-(b*r/4)^3))/((1-b*r/4)^3)-a*r^2;pp = 2*(p-r/3)/(G/3);psi=sqrt(sign(pp)*pp);fxs=fxs-G*w(k)*psi*cx(k);fys=fys-G*w(k)*psi*cy(k);end120D.3. Periodic boundary conditionsr=rho(i,j);p = (r*R*T*(1+(b*r/4)+(b*r/4)^2-(b*r/4)^3))/((1-b*r/4)^3)-a*r^2;pp = 2*(p-r/3)/(G/3);psi=sqrt(sign(pp)*pp);fx(i,j)=psi*fxs;fy(i,j)=psi*fys;endendendD.3 Periodic boundary conditionsPeriodic boundary conditions are implemented over each boundary lattice node, avoiding anyredundancies. This requires that only discrete velocities leaving the domain are considered forre-initialization on the opposite side of the domain, and that only one periodic boundary existsfor each cell.%Periodic boundary condition%left boundaryn(1,:,1)=n(nx,:,1);n(1,:,5)=n(nx,:,5);n(1,:,8)=n(nx,:,8);%right hand boundaryn(nx,:,3)=n(1,:,3);n(nx,:,7)=n(1,:,7);n(nx,:,6)=n(1,:,6);%bottom boundaryn(:,1,2)=n(:,ny,2);n(:,1,5)=n(:,ny,5);n(:,1,6)=n(:,ny,6);%Top boundaryfor i=2:nx-1n(i,ny,4)=n(i,1,4);n(i,ny,8)=n(i,1,8);n(i,ny,7)=n(i,1,7);endD.4 Collision operationThe single-relaxation collision operation for a Shan-Chen forcing scheme is representedin two-dimensions and three-dimensions.%% Collision operation% omega : represents 1/tau% ux: mean velocity in the x-direction% uy: mean velocity in the y-direction% vx: discrete velocity components (x-directed)% vy: discrete velocity components (y-directed)% w: lattice weight% forcx, forcy: x- and y-directed dimensionless external force density.function [n,uxf,uyf]=collision(nx,ny,ux,uy,vx,vy,omega,n,rho,w,forcx,forcy)neq = zeros(size(n));121D.4. Collision operationuxf = (1/omega)*forcx./rho+ux; % velocity with imparted momentum: x -->uyf = (1/omega)*forcy./rho+uy; % velocity with imparted momentum: y ^for j=1:nyfor i=1:nxt1=uxf(i,j)*uxf(i,j)+uyf(i,j)*uyf(i,j); % norm of velocity squaredfor m=1:9t2=uxf(i,j)*vx(m)+uyf(i,j)*vy(m);% equilibrium conditionneq(i,j,m)=rho(i,j)*w(m)*(1.0+3.0*t2+4.5*t2*t2-1.5*t1);% collision operationn(i,j,m)=(1.-omega)*n(i,j,m)+omega*neq(i,j,m);endendendendMulti-relaxation time collision operators require conversion between the moment spaceand the velocity space. Such conversion is accomplished using a transformation matrix.% Multi-relaxation-time collision functionfunction [n]=collision(nx,ny,ux,uy,n,rho,M,Minv,Deltam,forcx,forcy)I = eye(size(M));% Calculate the equilibrium condition in moment spacenmeq(:,:,1)= rho.*(-2+3*(ux.*ux+uy.*uy));nmeq(:,:,2)= rho.*( 1-3*(ux.*ux+uy.*uy));nmeq(:,:,3)= rho.* ux;nmeq(:,:,4)=-rho.* ux;nmeq(:,:,5)= rho.* uy;nmeq(:,:,6)=-rho.* uy;nmeq(:,:,7)= rho.*(ux.*ux-uy.*uy);nmeq(:,:,8)= rho.* ux.*uy;nmeq(:,:,9)= rho;%smeq=zeros(size(fmeq));smeq(:,:,1)=6*(ux.*forcx+uy.*forcy);smeq(:,:,2)=-6*(ux.*forcx+uy.*forcy);smeq(:,:,3)=forcx;smeq(:,:,4)=-forcx;smeq(:,:,5)=forcy;smeq(:,:,6)=-forcy;smeq(:,:,7)=2*(ux.*forcx-uy.*forcy);smeq(:,:,8)=ux.*forcy+uy.*forcx;smeq(:,:,9)=0;% transform ’n’ into moment space as ’nmom’for i=1:nxfor j=1:nyfor q=1:9smf=0.0;for m=1:9smf=smf+M(q,m)*n(i,j,m);endnmom(i,j,q)=smf;endendend% collision of ’nmom’ in momentum space with forcing term added.for j=1:ny122D.5. Propagation operationfor i=1:nxfor q=1:9ssmb=0.0;for m=1:9ssmb=ssmb-Deltam(q,m)*(nmom(i,j,m)-nmeq(i,j,m))+ ...(I(q,m)-(1/2)*Deltam(q,m))*smeq(i,j,m);% with force addedendnmom(i,j,q)=ssmb;endendend% transform ’nmom’ back into velocity space as ’n’for i=1:nxfor j=1:nyfor q=1:9smf=0.0;for m=1:9smf=smf+Minv(q,m)*nmom(i,j,m);endn(i,j,q)=n(i,j,q)+smf;endendendendD.5 Propagation operationThe propagation/streaming operation is identical in each method, and is greatly simplified bycircularly shifting of the discrete outcomes of nm(x, t). In this manner, the propagation/streamingoperation program can be generalized to represent any discrete velocity set.% Propagation/streaming:function [n]=stream(n,vx,vy)for m=1:8n(:,:,m)=circshift(squeeze(n(:,:,m)), [vx(m),vy(m)] );endend123
Thesis/Dissertation
eng
Electrical Engineering
Vancouver : University of British Columbia Library
Attribution-NonCommercial-NoDerivatives 4.0 International
Graduate
Chau, Kenneth
Simulations of optical radiation pressure on deformable media
Text
Still Image