Wave Optics Module: User's Guide
Wave Optics Module: User's Guide
Wave Optics Module: User's Guide
User’s Guide
Wave Optics Module User’s Guide
© 1998–2018 COMSOL
Protected by patents listed on www.comsol.com/patents, and U.S. Patents 7,519,518; 7,596,474;
7,623,991; 8,457,932; 8,954,302; 9,098,106; 9,146,652; 9,323,503; 9,372,673; and 9,454,625. Patents
pending.
This Documentation and the Programs described herein are furnished under the COMSOL Software License
Agreement (www.comsol.com/comsol-license-agreement) and may be used or copied only under the terms
of the license agreement.
COMSOL, the COMSOL logo, COMSOL Multiphysics, COMSOL Desktop, COMSOL Server, and
LiveLink are either registered trademarks or trademarks of COMSOL AB. All other trademarks are the
property of their respective owners, and COMSOL AB and its subsidiaries and products are not affiliated
with, endorsed by, sponsored by, or supported by those trademark owners. For a list of such trademark
owners, see www.comsol.com/trademarks.
Version: COMSOL 5.4
Contact Information
Visit the Contact COMSOL page at www.comsol.com/contact to submit general
inquiries, contact Technical Support, or search for an address and phone number. You can
also visit the Worldwide Sales Offices page at www.comsol.com/contact/offices for
address and contact information.
If you need to contact Support, an online request form is located at the COMSOL Access
page at www.comsol.com/support/case. Other useful links include:
Chapter 1: Introduction
Simplifying Geometries 25
2D Models . . . . . . . . . . . . . . . . . . . . . . . . . 25
3D Models . . . . . . . . . . . . . . . . . . . . . . . . . 26
Using Efficient Boundary Conditions . . . . . . . . . . . . . . . 27
Applying Electromagnetic Sources . . . . . . . . . . . . . . . . 28
Meshing and Solving . . . . . . . . . . . . . . . . . . . . . . 28
CONTENTS |3
Maxwell’s Equations 41
Introduction to Maxwell’s Equations . . . . . . . . . . . . . . . 41
Constitutive Relations . . . . . . . . . . . . . . . . . . . . . 42
Potentials. . . . . . . . . . . . . . . . . . . . . . . . . . 44
Electromagnetic Energy . . . . . . . . . . . . . . . . . . . . 44
Material Properties . . . . . . . . . . . . . . . . . . . . . . 45
About the Optical Materials Database . . . . . . . . . . . . . . . 47
Boundary and Interface Conditions . . . . . . . . . . . . . . . . 47
Phasors . . . . . . . . . . . . . . . . . . . . . . . . . . 48
Special Calculations 50
S-Parameter Calculations . . . . . . . . . . . . . . . . . . . . 50
Far-Field Calculations Theory . . . . . . . . . . . . . . . . . . 53
References . . . . . . . . . . . . . . . . . . . . . . . . . 54
Electromagnetic Quantities 64
4 | CONTENTS
Electromagnetic Waves, Frequency Domain Interface . . . . . . . 75
Wave Equation, Electric . . . . . . . . . . . . . . . . . . . . 77
Initial Values . . . . . . . . . . . . . . . . . . . . . . . . 82
External Current Density . . . . . . . . . . . . . . . . . . . 82
Far-Field Domain . . . . . . . . . . . . . . . . . . . . . . . 82
Far-Field Calculation . . . . . . . . . . . . . . . . . . . . . 83
Polarization . . . . . . . . . . . . . . . . . . . . . . . . . 84
Perfect Electric Conductor . . . . . . . . . . . . . . . . . . . 84
Perfect Magnetic Conductor . . . . . . . . . . . . . . . . . . 85
Port . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
Circular Port Reference Axis . . . . . . . . . . . . . . . . . . 93
Diffraction Order . . . . . . . . . . . . . . . . . . . . . . 93
Periodic Port Reference Point . . . . . . . . . . . . . . . . . . 95
Electric Field . . . . . . . . . . . . . . . . . . . . . . . . 96
Magnetic Field . . . . . . . . . . . . . . . . . . . . . . . . 97
Scattering Boundary Condition . . . . . . . . . . . . . . . . . 97
Reference Point . . . . . . . . . . . . . . . . . . . . . . 100
Impedance Boundary Condition . . . . . . . . . . . . . . . . 101
Surface Current Density . . . . . . . . . . . . . . . . . . . 103
Surface Magnetic Current Density . . . . . . . . . . . . . . . 103
Transition Boundary Condition . . . . . . . . . . . . . . . . 104
Periodic Condition . . . . . . . . . . . . . . . . . . . . . 106
Magnetic Current . . . . . . . . . . . . . . . . . . . . . 107
Edge Current . . . . . . . . . . . . . . . . . . . . . . . 108
Electric Point Dipole . . . . . . . . . . . . . . . . . . . . 108
Magnetic Point Dipole . . . . . . . . . . . . . . . . . . . . 108
Line Current (Out-of-Plane) . . . . . . . . . . . . . . . . . 109
CONTENTS |5
Wave Equations . . . . . . . . . . . . . . . . . . . . . . 120
Initial Values . . . . . . . . . . . . . . . . . . . . . . . 122
Electric Current Density . . . . . . . . . . . . . . . . . . . 123
Magnetic Current Density . . . . . . . . . . . . . . . . . . 123
Electric Field . . . . . . . . . . . . . . . . . . . . . . . 123
Perfect Electric Conductor . . . . . . . . . . . . . . . . . . 124
Magnetic Field . . . . . . . . . . . . . . . . . . . . . . . 124
Perfect Magnetic Conductor . . . . . . . . . . . . . . . . . 124
Surface Current Density . . . . . . . . . . . . . . . . . . . 125
Scattering Boundary Condition . . . . . . . . . . . . . . . . 125
Flux/Source . . . . . . . . . . . . . . . . . . . . . . . . 126
Background Field . . . . . . . . . . . . . . . . . . . . . . 126
Far-Field Domain . . . . . . . . . . . . . . . . . . . . . . 127
Far-Field Calculation . . . . . . . . . . . . . . . . . . . . 127
6 | CONTENTS
Interface 156
The Equations . . . . . . . . . . . . . . . . . . . . . . . 156
In-plane E Field or In-plane H Field . . . . . . . . . . . . . . . 160
Fluxes as Dirichlet Boundary Conditions . . . . . . . . . . . . . 161
Absorbing Layers . . . . . . . . . . . . . . . . . . . . . . 162
Chapter 5: Glossary
CONTENTS |7
8 | CONTENTS
1
Introduction
This guide describes the Wave Optics Module, an optional add-on package for
COMSOL Multiphysics® designed to assist you to set up and solve electromagnetic
wave problems at optical frequencies.
This chapter introduces you to the capabilities of this module. A summary of the
physics interfaces and where you can find documentation and model examples is
also included. The last section is a brief overview with links to each chapter in this
guide.
In this chapter:
9
About the Wave Optics Module
These topics are included in this section:
The Wave Optics Module solves problems in the field of electromagnetic waves at
optical frequencies (corresponding to wavelengths in the nano- to micrometer range).
The underlying equations for electromagnetics are automatically available in all of the
physics interfaces — a feature unique to COMSOL Multiphysics. This also makes
nonstandard modeling easily accessible.
The module is useful for simulations and design of optical applications in virtually all
areas where you find electromagnetic waves, such as:
• Optical fibers
• Photonic waveguides
• Photonic crystals
• Nonlinear optics
10 | CHAPTER 1: INTRODUCTION
• Laser resonator design
• Active devices in photonics
The physics interfaces cover the following types of electromagnetics field simulations
and handle time-harmonic, time-dependent, and eigenfrequency/eigenmode
problems:
Material properties include inhomogeneous and fully anisotropic materials, media with
gains or losses, and complex-valued material properties. In addition to the standard
postprocessing features, the module supports direct computation of S-parameters and
far-field radiation patterns. You can add ports with a wave excitation with specified
power level and mode type, and add PMLs (perfectly matched layers) to simulate
electromagnetic waves that propagate into an unbounded domain. For time-harmonic
simulations, you can use the scattered wave or the total wave.
Both the RF and the Wave Optics Module can handle high-frequency electromagnetic
wave simulations. However, with the Wave Optics Module you can do time-harmonic
simulations of domains that are much larger than the wavelength. This situation is
typical for optical phenomena, components, and systems. Due to the relatively weak
coupling between waves in optical materials, the interaction lengths are often much
larger than the wavelength. This applies to linear couplers, like directional couplers and
fiber Bragg gratings, and nonlinear phenomena, like second harmonic generation, self-
phase modulation, and so forth. With the Wave Optics Module, these kinds of
problems are directly addressable, without huge computer memory requirements.
Heat Transfer
Electromagnetic Heating
12 | CHAPTER 1: INTRODUCTION
PHYSICS INTERFACE ICON TAG SPACE AVAILABLE STUDY TYPE
DIMENSION
Optics
Wave Optics
However, if the Frequency Domain or Wavelength Domain studies, available with the
Electromagnetic Waves, Frequency Domain and the Electromagnetic Waves, Beam
Envelopes interfaces, are used, this allows you to efficiently simplify and assume that
all variations in time occur as sinusoidal signals. Then the problem is time harmonic
and in the frequency domain. Thus you can formulate it as a stationary problem with
complex-valued solutions. The complex value represents both the amplitude and the
phase of the field, while the frequency is specified as a scalar model input, usually
provided by the solver. This approach is useful because, combined with Fourier
analysis, it applies to all periodic signals with the exception of nonlinear problems.
14 | CHAPTER 1: INTRODUCTION
Examples of typical frequency or wavelength domain simulations are wave-
propagation problems.
For nonlinear problems you can apply a Frequency Domain or Wavelength Domain
study after a linearization of the problem, which assumes that the distortion of the
sinusoidal signal is small. You can also couple waves at different frequencies, for
example in applications like second harmonic generation, by coupling several physics
interfaces, defined for the different frequencies, using weak expression coupling terms.
Use a Time Dependent study when the nonlinear influence is strong, or if you are
interested in the harmonic distortion of a sine signal. It can also be more efficient to
use a time dependent study if you have a periodic input with many harmonics, like a
square-shaped signal.
For many optical applications the propagation length is much longer than the
wavelength. For instance, a typical optical wavelength is 1 μm, but the propagation
length can easily be on the millimeter to centimeter scale. To apply the Frequency
Domain interface to this kind of problems, requires a large amount of available
memory. However, many problems are such that the electric field can be factored into
a slowly varying amplitude factor and a rapidly varying phase factor. The
Electromagnetic Waves, Beam Envelopes interface is based on this assumption. Thus,
this physics interface assumes a prescribed rapidly varying phase factor and solves for
the slowly varying amplitude factor. Thereby it can be used for solving problems
extending over domains that are a large number of wavelengths long, without
requiring the use of large amounts of memory.
Even after a model is defined, you can edit input data, equations, boundary conditions,
geometry — the equations and boundary conditions are still available through
associative geometry — and mesh settings. You can restart the solver, for example,
using the existing solution as the initial condition or initial guess. It is also easy to add
another physics interface to account for a phenomenon not previously described in a
model.
16 | CHAPTER 1: INTRODUCTION
THE DOCUMENTATION AND ONLINE HELP
The COMSOL Multiphysics Reference Manual describes the core physics interfaces
and functionality included with the COMSOL Multiphysics license. This book also has
instructions about how to use COMSOL Multiphysics and how to access the
electronic Documentation and Help content.
• In the Model Builder or Physics Builder click a node or window and then
press F1.
• On the main toolbar, click the Help ( ) button.
• From the main menu, select Help>Help.
• Press Ctrl+F1.
• From the File menu select Help>Documentation ( ).
• Press Ctrl+F1.
• On the main toolbar, click the Documentation ( ) button.
• From the main menu, select Help>Documentation.
Once the Application Libraries window is opened, you can search by name or browse
under a module folder name. Click to view a summary of the model or application and
its properties, including options to open it or its associated PDF document.
To include the latest versions of model examples, from the Help menu
select ( ) Update COMSOL Application Library.
18 | CHAPTER 1: INTRODUCTION
CONTACTING COMSOL BY EMAIL
For general product information, contact COMSOL at [email protected].
The chapter also contains a review of the basic theory of electromagnetics, starting
with Maxwell’s Equations, and the theory for some Special Calculations: S-parameters,
and far-field analysis. There is also a list of Electromagnetic Quantities with the SI units
and symbols.
OPTICS
Wave Optics Interfaces chapter describes:
20 | CHAPTER 1: INTRODUCTION
• The Electromagnetic Waves, Time Explicit Interface, which solves a transient wave
equation for both the electric and magnetic fields.
• The Electromagnetic Waves, Beam Envelopes Interface, which analyzes frequency
domain electromagnetic waves, and uses time-harmonic and eigenfrequency
studies, boundary mode analysis, and frequency domain, modal studies.
HEAT TRANSFER
The Laser Heating Interface is used to model electromagnetic heating for systems and
devices where the electric field amplitude varies slowly on a wavelength scale. This
multiphysics interface adds an Electromagnetic Waves, Beam Envelopes interface and
a Heat Transfer in Solids interface.
The goal of this chapter is to familiarize you with the modeling procedure in the
Wave Optics Module. A number of models available in the Applications Libraries
also illustrate the different aspects of the simulation process.
In this chapter:
23
Preparing for Wave Optics Modeling
Several modeling topics are described in this section that might not be found in
ordinary textbooks on electromagnetic theory.
Increasing the complexity of a model to make it more accurate usually makes it more
expensive to simulate. A complex model is also more difficult to manage and interpret
than a simple one. Keep in mind that it can be more accurate and efficient to use several
simple models instead of a single, complex one.
In this section:
• 2D Models
• 3D Models
• Using Efficient Boundary Conditions
• Applying Electromagnetic Sources
• Meshing and Solving
2D Models
The text below is a guide to some of the common approximations made for 2D
models. Remember that the modeling in 2D usually represents some 3D geometry
under the assumption that nothing changes in the third dimension or that the field has
a prescribed propagation component in the third dimension.
CARTESIAN COORDINATES
In this case a cross section is viewed in the xy-plane of the actual 3D geometry. The
geometry is mathematically extended to infinity in both directions along the z-axis,
assuming no variation along that axis or that the field has a prescribed wave vector
component along that axis. All the total flows in and out of boundaries are per unit
length along the z-axis. A simplified way of looking at this is to assume that the
geometry is extruded one unit length from the cross section along the z-axis. The total
flow out of each boundary is then from the face created by the extruded boundary (a
boundary in 2D is a line).
There are usually two approaches that lead to a 2D cross-section view of a problem.
The first approach is when it is known that there is no variation of the solution in one
SIMPLIFYING GEOMETRIES | 25
particular dimension. The second approach is when there is a problem where the
influence of the finite extension in the third dimension can be neglected.
When using the axisymmetric versions, the horizontal axis represents the
radial (r) direction and the vertical axis the z direction, and the geometry
in the right half-plane (that is, for positive r only) must be created.
POLARIZATION IN 2D
In addition to selecting 2D or 2D axisymmetry when you start building the model, the
physics interfaces (The Electromagnetic Waves, Frequency Domain Interface, The
Electromagnetic Waves, Transient Interface, or The Electromagnetic Waves, Beam
Envelopes Interface) in the Model Builder offers a choice in the Components settings
section. The available choices are Out-of-plane vector, In-plane vector, and
Three-component vector. This choice determines what polarizations can be handled.
For example, as you are solving for the electric field, a 2D TM (out-of-plane H field)
model requires choosing In-plane vector as then the electric field components are in
the modeling plane.
3D Models
Although COMSOL Multiphysics fully supports arbitrary 3D geometries, it is
important to simplify the problem. This is because 3D models often require more
computer power, memory, and time to solve. The extra time spent on simplifying a
model is probably well spent when solving it. Below are a few issues that need to be
addressed before starting to implement a 3D model in this module.
• Check if it is possible to solve the problem in 2D. Given that the necessary
approximations are small, the solution is more accurate in 2D, because a much
denser mesh can be used.
• Many models extend to infinity or can have regions where the solution only
undergoes small changes. This problem is addressed in two related steps. First, the
geometry needs to be truncated in a suitable position. Second, a suitable boundary
condition needs to be applied there. For static and quasi-static models, it is often
possible to assume zero fields at the open boundary, provided that this is at a
sufficient distance away from the sources. For radiation problems, special
low-reflecting boundary conditions need to be applied. This boundary should be in
the order of a few wavelengths away from any source.
A more accurate option is to use perfectly matched layers (PMLs). PMLs are layers
that absorbs all radiated waves with small reflections.
• Replace thin layers with boundary conditions where possible. There are several types
of boundary conditions in COMSOL Multiphysics suitable for such replacements.
SIMPLIFYING GEOMETRIES | 27
For example, replace materials with high conductivity by the perfect electric
conductor (PEC) boundary condition.
• Use boundary conditions for known solutions. For example, an antenna aperture
can be modeled as an equivalent surface current density on a 2D face (boundary) in
a 3D model.
• The first is the variation in the solution due to geometrical factors. The mesh
generator automatically generates a finer mesh where there is a lot of fine
geometrical details. Try to remove such details if they do not influence the solution,
because they produce a lot of unnecessary mesh elements.
• The second is the skin effect or the field variation due to losses. It is easy to estimate
the skin depth from the conductivity, permeability, and frequency. At least two linear
elements per skin depth are required to capture the variation of the fields. If the skin
depth is not studied or a very accurate measure of the dissipation loss profile is not
needed, replace regions with a small skin depth with a boundary condition, thereby
SOLVERS
In most cases the solver sequence generated by COMSOL Multiphysics can be used.
The choice of solver is optimized for the typical case for each physics interface and
study type in this module. However, in special cases tuning the solver settings can be
required. This is especially important for 3D problems because they can require a large
amount of memory.
• Meshing
• Studies and Solvers
SIMPLIFYING GEOMETRIES | 29
Periodic Boundary Conditions
The Wave Optics Module has a dedicated Periodic Condition. The periodic condition
can identify simple mappings on plane source and destination boundaries of equal
shape. The destination can also be rotated with respect to the source. There are three
types of periodic conditions available (only the first two for transient analysis):
• Continuity — The tangential components of the solution variables are equal on the
source and destination.
• Antiperiodicity — The tangential components have opposite signs.
• Floquet periodicity — There is a phase shift between the tangential components.
The phase shift is determined by a wave vector and the distance between the source
and destination. Floquet periodicity is typically used for models involving plane
waves interacting with periodic structures.
Periodic boundary conditions must have compatible meshes. This can be done
automatically by enabling the Physics-control mesh in the setting for The
Electromagnetic Waves, Frequency Domain Interface or by manually setting up the
correct mesh sequence
To learn how to use the Copy Mesh feature to ensure that the mesh on
the destination boundary is identical to that on the source boundary, see
Band-Gap Analysis of a Photonic Crystal: Application Library path
Wave_Optics_Module/Gratings_and_Metamaterials/
bandgap_photonic_crystal.
∇ × ( μ –1 ∇ × E ) – ω 2 εc E = 0
The known background field becomes a source term and the scattered field
formulation thus solves for the relative electric field. A linearly polarized plane wave
background field, a paraxial-approximate Gaussian beam, or a user-defined
background field can be specified. When solving the scattered field formulation the
total, the background, and the relative electric fields are available. The relative field is
the difference between the background field and the total field. It is the relative field
that contributes to the far-field calculation. For more information about the Far-Field
computation, see Far-Field Calculations Theory. The benefit to this approach is that if
the background field is much larger in magnitude than the scattered field, the accuracy
of the simulation improves if the relative field is solved for. Another advantage is that
is becomes very easy to set up a perfectly matched layer surrounding the homogeneous
medium modeling domain.
The drawback to the this approach is that the relative field requires some careful
interpretation. The relative electric field can conceptually be decomposed into:
The Escattered component is the scattered field from object. This is the field that is of
interest in a scattering problem. However, the relative field may also consist of a
component that represents a correction to the background field and a cancellation of
the background field. The Ecorrection component can be nonzero when the
background field does not exactly satisfy Maxwell’s equations, such as when the
paraxial Gaussian beam approximation is used for a tightly focused beam. For more
information about the Gaussian beam theory, see Gaussian Beams as Background
Fields. The Ecancellation component will be nonzero and equal to −Ebackground
wherever the total field should be zero, such as in the interior of any perfectly shielded
An alternative of using the scattered-field formulation, is to use ports with the Activate
slit condition on interior port setting enabled. Then the domain can be excited by the
interior port and the outgoing field can be absorbed by perfectly matched layers. For
more information about the Port feature and the Activate slit condition on interior port
setting, see Port Properties.
In this section:
The Far-Field Domain and the Far-Field Calculation nodes get their selections
automatically, if the Perfectly Matched Layer (PML) feature has been defined before
adding the Far-Field Domain feature.
For each variable name entered, the software generates functions and variables, which
represent the vector components of the far electric field. The names of these variables
For example, the name Efar is entered and the geometry is Cartesian with the
independent variables x, y, and z, the generated variables get the names Efarx, Efary,
and Efarz.
If, on the other hand, the geometry is axisymmetric with the independent variables r,
phi, and z, the generated variables get the names Efarr, Efarphi, and Efarz.
In 2D, the software only generates the variables for the nonzero field components. The
physics interface name also appears in front of the variable names so they can vary, but
typically look something like ewfd.Efarz and so forth.
To each of the generated variables, there is a corresponding function with the same
name. This function takes the vector components of the evaluated far-field direction as
arguments.
The vector components also can be interpreted as a position. For example, assume that
the variables dx, dy, and dz represent the direction in which the far electric field is
evaluated.
The expression
Efarx(dx,dy,dz)
gives the value of the far electric field in this direction. To give the direction as an angle,
use the expression
Efarx(sin(theta)*cos(phi),sin(theta)*sin(phi),cos(theta))
where the variables theta and phi are defined to represent the angular direction
( θ, ϕ ) in radians. The magnitude of the far field and its value in dB are also generated
as the variables normEfar and normdBEfar, respectively.
The main advantage with the Radiation Pattern plot, as compared to making a Line
Graph, is that the unit circle/sphere that you use for defining the plot directions, is not
part of your geometry for the solution. Thus, the number of plotting directions is
decoupled from the discretization of the solution domain.
be spherical for 3D and circular for 2D axisymmetric components and its center to be
on the origin.
23D far-field norm functions in 2D axisymmetric geometry are available in these cases:
• Antenna models using circular port excitation with a positive azimuthal mode
number
• Scattered field analysis excited by the predefined circularly polarized plane wave type
The function can be used in a 3D Radiation Pattern plot, where the input argument of
the function must be same as the Azimuth angle variable in the Evaluation section in the
settings window.
The suffix of a function name varies based on the circular port mode type, port mode
number and azimuthal mode number in the physics interface. For example, when
using azimuthal mode number 1 in the physics interface and transverse electric (TE)
When the function is used in a radiation pattern plot under a 1D or a polar plot group,
the value of input argument defines the plotting plane regardless of the normal and
reference direction in the evaluation section in the settings window. For example,
norm3DEfar_TE12(0)evaluates the norm of the electric far field for the TE12 mode
for 0 degree azimuthal angle. This is equivalent to plotting the this variable on the
xz-plane. Similarly, norm3DEfar_TE12(pi/2) is the evaluation at 90 degree azimuthal
angle, which is equivalent to plotting the variable on the yz-plane.
The 3D far-field norm, the linear superposition of the positive and negative azimuthal
modes scaled by 0.5, is
2 2 2
E r cos mφ + E φ sin m φ + E z cos mφ ,
nz
sin ----- ( 2πd z sin θ cos φ + α z )
2
× ----------------------------------------------------------------------------- ,
2πd sin θ cos φ + α
sin ---------------------------------------------------
z z
2
The uniform two dimensional array factor is simpler than the three dimensional
version, as the third, the z-component factor, is unity.
4In the postprocessing expression context menu, far-field functions are available under
Gain
The antenna realized gain is defined as
2
4πU normEfar
G real ized = ------------ = ----------------------------------
P in 60P in
∗ 2
where U is the radiation intensity, Re ( E far × H far ) ⁄ 2 = normEfar ⁄ 240π and
Pin is the total input power.
2
normEfar
G = ----------------------------------
60P delivered
2
where the delivered power, Pdelivered is P in ( 1 – S 11 ) . The gain is available only
when the S-parameter calculation is valid, that is, for the single port excitation case.
The equations can be formulated in differential or integral form. The differential form
are presented here, because it leads to differential equations that the finite element
method can handle. For general time-varying fields, Maxwell’s equations can be
written as
MAXWELL’S EQUATIONS | 41
∂D
∇ × H = J + -------
∂t
∂------
B-
∇×E = –
∂t
∇⋅D = ρ
∇⋅B = 0
The first two equations are also referred to as Maxwell-Ampère’s law and Faraday’s
law, respectively. Equation three and four are two forms of Gauss’ law, the electric and
magnetic form, respectively.
∂ρ
∇⋅J = – ------
∂t
Out of the five equations mentioned, only three are independent. The first two
combined with either the electric form of Gauss’ law or the equation of continuity
form such an independent system.
Constitutive Relations
To obtain a closed system, the constitutive relations describing the macroscopic
properties of the medium, are included. They are given as
D = ε0 E + P
B = μ0 ( H + M )
J = σE
Here ε0 is the permittivity of vacuum, μ0 is the permeability of vacuum, and σ the
electrical conductivity. In the SI system, the permeability of a vacuum is chosen to be
4π·10−7 H/m. The velocity of an electromagnetic wave in a vacuum is given as c0 and
the permittivity of a vacuum is derived from the relation
1 – 12 1 –9
ε 0 = ----------
2
= 8.854 ⋅ 10 F/m ≈ --------- ⋅ 10 F/m
c0 μ0 36π
The electric polarization vector P describes how the material is polarized when an
electric field E is present. It can be interpreted as the volume density of electric dipole
moments. P is generally a function of E. Some materials can have a nonzero P also
when there is no electric field present.
For linear materials, the polarization is directly proportional to the electric field,
P = ε0χeE, where χe is the electric susceptibility. Similarly in linear materials, the
magnetization is directly proportional to the magnetic field, M = χmH, where χm is the
magnetic susceptibility. For such materials, the constitutive relations can be written
D = ε 0 ( 1 + χ e )E = ε 0 ε r E = εE
B = μ 0 ( 1 + χ m )H = μ 0 μ r H = μH
The parameter εr is the relative permittivity and μr is the relative permeability of the
material. These are usually scalar properties but they can, for a general anisotropic
material, be 3-by-3 tensors. The properties ε and μ (without subscripts) are the
permittivity and permeability of the material.
D = ε0 εr E + Dr
The field Dr is the remanent displacement, which is the displacement when no electric
field is present.
Similarly, a generalized form of the constitutive relation for the magnetic field is
B = μ0 μr H + Br
where Br is the remanent magnetic flux density, which is the magnetic flux density
when no magnetic field is present.
e
J = σE + J
MAXWELL’S EQUATIONS | 43
Potentials
Under certain circumstances it can be helpful to formulate the problems in terms of
the electric scalar potential V and the magnetic vector potential A. They are given by
the equalities
B = ∇×A
∂A
E = – ∇V – -------
∂t
The defining equation for the magnetic vector potential is a direct consequence of the
magnetic Gauss’ law. The electric potential results from Faraday’s law.
Electromagnetic Energy
The electric and magnetic energies are defined as
D T
∂D
We = V 0 E ⋅ dD dV = V 0 E ⋅ -------
∂t
dt dV
B T
∂B
Wm = 0 H ⋅ dB dV = - dt dV
0 H ⋅ ------
V V ∂t
The time derivatives of these expressions are the electric and magnetic power
∂D
Pe = V E ⋅ -------
∂t
dV
∂B
Pm = V H ⋅ ------
∂t
- dV
These quantities are related to the resistive and radiative energy, or energy loss,
through Poynting’s theorem (Ref. 3)
∂D ∂B
– V E ⋅ -------
∂t
+ H ⋅ ------- dV = J ⋅ E dV + ( E × H ) ⋅ n dS
∂t V °S
where V is the computation domain and S is the closed boundary of V.
The first term on the right-hand side represents the resistive losses,
Ph = V J ⋅ E dV
The second term on the right-hand side of Poynting’s theorem represents the radiative
losses,
Pr =
°S ( E × H ) ⋅ n dS
The quantity S = E × H is called the Poynting vector.
Under the assumption the material is linear and isotropic, it holds that
∂D ∂E
E ⋅ ------- = εE ⋅ ------- = ∂ --- εE ⋅ E
1
∂t ∂t ∂t 2
∂B ∂B
H ⋅ ------- = --- B ⋅ ------- = ∂ ------- B ⋅ B
1 1
∂t μ ∂t ∂ t 2μ
By interchanging the order of differentiation and integration (justified by the fact that
the volume is constant and the assumption that the fields are continuous in time), this
equation results:
∂
V --2- εE ⋅ E + ------
- B ⋅ B dV =
1 1
–
∂t 2μ V J ⋅ E dV + °S ( E × H ) ⋅ n dS
The integrand of the left-hand side is the total electromagnetic energy density
1 1
w = w e + w m = --- εE ⋅ E + ------- B ⋅ B
2 2μ
Material Properties
Until now, there has only been a formal introduction of the constitutive relations.
These seemingly simple relations can be quite complicated at times. There are four
main groups of materials where they require some considerations. A given material can
belong to one or more of these groups. The groups are:
• Inhomogeneous materials
• Anisotropic materials
• Nonlinear materials
• Dispersive materials
MAXWELL’S EQUATIONS | 45
The least complicated of the groups above is that of the inhomogeneous materials. An
inhomogeneous medium is one where the constitutive parameters vary with the space
coordinates, so that different field properties prevail at different parts of the material
structure.
For anisotropic materials, the field relations at any point are different for different
directions of propagation. This means that a 3-by-3 tensor is required to properly
define the constitutive relations. If this tensor is symmetric, the material is often
referred to as reciprocal. In these cases, the coordinate system can be rotated in such
a way that a diagonal matrix is obtained. If two of the diagonal entries are equal, the
material is uniaxially anisotropic. If none of the elements have the same value, the
material is biaxially anisotropic (Ref. 2). An example where anisotropic parameters
are used is for the permittivity in crystals (Ref. 2).
Finally, dispersion describes changes in the velocity of the wave with wavelength. In
the frequency domain, dispersion is expressed by a frequency dependence in the
constitutive laws.
The physics-specific domain material properties are by default taken from the material
specification. The material properties are inputs to material laws or constitutive
relations that are defined on the feature level below the physics interface node in the
model tree. There is one editable default domain feature (wave equation) that initially
represents a linear isotropic material. Domains with different material laws are specified
by adding additional features. Some of the domain parameters can either be a scalar or
a matrix (tensor) depending on whether the material is isotropic or anisotropic.
The real and the imaginary parts of the refractive index are available as interpolation
functions for a number of organic and inorganic materials, including many glasses and
some semiconductors.
For an example of how to use the Optical materials database, see Optical
Scattering Off a Gold Nanosphere: Application Library path
Wave_Optics_Module/Optical_Scattering/scattering_nanosphere
n2 × ( E1 – E2 ) = 0
n2 ⋅ ( D1 – D2 ) = ρs
n2 × ( H1 – H2 ) = Js
n2 ⋅ ( B1 – B2 ) = 0
where ρs and Js denote surface charge density and surface current density,
respectively, and n2 is the outward normal from medium 2. Of these four conditions,
only two are independent. One of the first and the fourth equations, together with one
of the second and third equations, form a set of two independent conditions.
A consequence of the above is the interface condition for the current density,
MAXWELL’S EQUATIONS | 47
∂ρ s
n 2 ⋅ ( J 1 – J 2 ) = – --------
∂t
–n2 × E2 = 0
–n2 × H2 = Js
–n2 ⋅ D2 = ρs
–n2 ⋅ B2 = 0
Phasors
Whenever a problem is time-harmonic the fields can be written in the form
ˆ
E ( r, t ) = E ( r ) cos ( ωt + φ )
Instead of using a cosine function for the time dependence, it is more convenient to
use an exponential function, by writing the field as
ˆ ˆ jφ jωt ˜ jωt
E ( r, t ) = E ( r ) cos ( ωt + φ ) = Re ( E ( r )e e ) = Re ( E ( r )e )
˜
The field E ( r ) is a phasor (phase vector), which contains amplitude and phase
information of the field but is independent of t. One thing that makes the use of
phasors suitable is that a time derivative corresponds to a multiplication by jω,
∂------
E- ˜ jωt
= Re ( jωE ( r )e )
∂t
This means that an equation for the phasor can be derived from a time-dependent
equation by replacing the time derivatives by a factor jω. All time-harmonic equations
˜
Re ( E ( r ) )
MAXWELL’S EQUATIONS | 49
Special Calculations
In this section:
• S-Parameter Calculations
• Far-Field Calculations Theory
• References
S-Parameter Calculations
For high-frequency problems, voltage is not a well-defined entity, and it is necessary
to define the scattering parameters (S-parameter) in terms of the electric field. To
convert an electric field pattern on a port to a scalar complex number corresponding
to the voltage in transmission line theory an eigenmode expansion of the
electromagnetic fields on the ports needs to be performed. Assume that an eigenmode
analysis has been performed on the ports 1, 2, 3, … and that the electric field patterns
E1, E2, E3, … of the fundamental modes on these ports are known. Further, assume
that the fields are normalized with respect to the integral of the power flow across each
port cross section, respectively. This normalization is frequency dependent unless
TEM modes are being dealt with. The port excitation is applied using the fundamental
eigenmode, the mode with subscript 1. The computed electric field Ec on the port
consists of the excitation plus the reflected field. That is, on the port boundary where
there is an incident wave, the computed field can be expanded in terms of the mode
fields as
Ec = E1 + Si1 Ei ,
i=1
Ec = Si1 Ei
i=1
The S-parameter for the mode with index k is then given by multiplying with the
conjugate of the mode field for mode k and integrating over the port boundary. Since
the mode fields for the different modes are orthogonal, the following relations are
obtained for the S-parameters
*
( E c ⋅ E 2 ) dA 2
port 2
S 21 = ----------------------------------------------
-
*
( E 2 ⋅ E 2 ) dA 2
port 2
*
( E c ⋅ E 3 ) dA 3
port 3
S 31 = ----------------------------------------------
-
*
( E 3 ⋅ E 3 ) dA 3
port 3
and so on. To get S22 and S12, excite port number 2 in the same way.
1 *
S av = --- Re ( E × H )
2
The amount of power flowing out of a port is given by the normal component of the
Poynting vector,
1 *
n ⋅ S av = n ⋅ --- Re ( E × H )
2
Below the cutoff frequency the power flow is zero, which implies that it is not possible
to normalize the field with respect to the power flow below the cutoff frequency. But
in this region the S-parameters are trivial and do not need to be calculated.
In the following subsections the power flow is expressed directly in terms of the electric
field for TE, TM, and TEM waves.
TE Waves
For TE waves it holds that
E = – Z TE ( n × H )
SPECIAL CALCULATIONS | 51
where ZTE is the wave impedance
ωμ
Z TE = -------
β
ω is the angular frequency of the wave, μ the permeability, and β the propagation
constant. The power flow then becomes
1 * 1 * 1 2
n ⋅ S av = --- n ⋅ Re ( E × H ) = – --- Re ( E ⋅ ( n × H ) ) = -------------- E
2 2 2Z TE
TM Waves
For TM waves it holds that
1
H = ----------- ( n × E )
Z TM
β
Z TM = -------
ωε
1 * 1 *
n ⋅ S av = --- n ⋅ Re ( E × H ) = --------------- ( n ⋅ Re ( E × ( n × E ) ) )
2 2Z TM
1 2
= --------------- n × E
2Z TM
TEM Waves
For TEM waves it holds that
1
H = --------------- ( n × E )
Z TEM
μ
Z TEM = ---
ε
1 * 1 2 1 2
n ⋅ S av = --- n ⋅ Re ( E × H ) = ------------------ n × E = ------------------ E
2 2Z TEM 2Z TEM
jk
E p = ------ r 0 × [ n × E – ηr 0 × ( n × H ) ] exp ( jkr ⋅ r 0 ) dS
4π
jk
Ep =
4π
λ ------ r 0 × [ n × E – ηr 0 × ( n × H ) ] exp ( jkr ⋅ r 0 ) dS
In both cases the integration is performed on a closed boundary. In the scattered field
formulation, where the total electric field is the sum of the background field and the
scattered field, the far-field only gets contributions from the scattered field, since the
contributions from the background field cancel out when integrated over all parts of
the closed boundary.
For scattering problems, the far field in COMSOL Multiphysics is identical to what in
physics is known as the “scattering amplitude”.
The antenna is located in the vicinity of the origin, while the far-field point p is taken
at infinity but with a well-defined angular position ( θ, ϕ ) .
• E and H are the fields on the “aperture”—the surface S enclosing the antenna.
• r0 is the unit vector pointing from the origin to the field point p. If the field points
lie on a spherical surface S', r0 is the unit normal to S'.
• n is the unit normal to the surface S.
• η is the impedance:
η = μ⁄ε
SPECIAL CALCULATIONS | 53
Thus the unit vector r0 can be interpreted as the direction defined by the angular
position ( θ, ϕ ) and Ep is the far field in this direction.
Because the far field is calculated in free space, the magnetic field at the far-field point
is given by
r0 × Ep
H p = -------------------
η0
The Poynting vector gives the power flow of the far field:
* 2
r 0 ⋅ S = r 0 ⋅ Re ( E p × H p ) ∼ E p
References
1. D.K. Cheng, Field and Wave Electromagnetics, 2nd ed., Addison-Wesley, 1991.
4. R.K. Wangsness, Electromagnetic Fields, 2nd ed., John Wiley & Sons, 1986.
S 11 S 12 . . S 1n
S 21 S 22 . . .
S = . . . . .
. . . . .
S n1 . . . S nn
where S11 is the voltage reflection coefficient at port 1, S21 is the voltage transmission
coefficient from port 1 to port 2, and so on. The time average power reflection/
transmission coefficients are obtained as |Sij |2.
S-Parameter Variables
This module automatically generates variables for the S-parameters. The port names
(use numbers for sweeps to work correctly) determine the variable names. If, for
example, there are two ports with the numbers 1 and 2 and Port 1 is the inport, the
software generates the variables S11 and S21. S11 is the S-parameter for the reflected
wave and S21 is the S-parameter for the transmitted wave. For convenience, two
variables for the S-parameters on a dB scale, S11dB and S21dB, are also defined using
the following relation:
S 11dB = 20 log 10 ( S 11 )
The model and physics interface names also appear in front of the variable names so
they can vary. The S-parameter variables are added to the predefined quantities in
appropriate plot lists.
In this section:
• Eigenfrequency Analysis
• Mode Analysis
Eigenfrequency Analysis
The eigenfrequency analysis solves for the eigenfrequency of a model. The
time-harmonic representation of the fields is more general and includes a complex
parameter in the phase
˜ jωt ˜ –λ t
E ( r, t ) = Re ( E ( r T )e ) = Re ( E ( r )e )
where the eigenvalue, (−λ) = −δ + jω, has an imaginary part representing the
eigenfrequency, and a real part responsible for the damping. It is often more common
to use the quality factor or Q-factor, which is derived from the eigenfrequency and
damping
ω
Q fact = ---------
2δ
These situations may require special treatment, especially since it can lead to “singular
matrix” or “undefined value” messages if not treated correctly. Under normal
circumstances, the automatically generated solver settings should handle the cases
described in the table above. However, the following discussion provide some
background to the problem of defining the eigenvalue linearization point. The
complication is not only the nonlinearity itself, it is also the way it enters the equations.
For example the impedance boundary conditions with nonzero boundary conductivity
has the term
ε 0 μ 0 μ rbnd
– ( – λ ) ------------------------------------------ ( n × ( n × H ) )
σ bnd
ε rbnd + -----------------
( – λ )ε 0
where (−λ) = −δ + jω. When the solver starts to solve the eigenfrequency problem it
linearizes the entire formulation with respect to the eigenvalue around a certain
linearization point. By default this linearization point is set to the value provided to the
Search for eigenvalues around field, for the three cases listed in the table above.
Normally, this should be a good value for the linearization point. For instance, for the
impedance boundary condition, this avoids setting the eigenvalue λ to zero in the
denominator in the equation above. For other cases than those listed in the table
above, the default linearization point is zero.
If the default values for the linearization point is not suitable for your particular
problem, you can manually provide a “good” linearization point for the eigenvalue
solver. Do this in the Eigenvalue node (not the Eigenfrequency node) under the Solver
Sequence node in the Study branch of the Model Builder. A solver sequence can be
generated first. In the Linearization Point section, select the Transform point check box
In many cases it is enough to specify a good linearization point and then solve the
problem once. If a more accurate eigenvalue is needed, an iterative scheme is necessary:
1 Specify that the eigenvalue solver only searches for one eigenvalue. Do this either
for an existing solver sequence in the Eigenvalue node or, before generating a solver
sequence, in the Eigenfrequency node.
2 Solve the problem with a “good” linearization point. As the eigenvalue shifts, use
the same value with the real part removed from the eigenvalue or, equivalently, use
the real part of the eigenfrequency.
3 Extract the eigenvalue from the solution and update the linearization point and the
shift.
4 Repeat until the eigenvalue does not change more than a desired tolerance.
• For a list of the studies available by physics interface, see The Wave
Optics Module Physics Interface Guide
• Studies and Solvers in the COMSOL Multiphysics Reference Manual
Mode Analysis
In mode analysis and boundary mode analysis COMSOL Multiphysics solves for the
propagation constant. The time-harmonic representation is almost the same as for the
eigenfrequency analysis, but with a known propagation in the out-of-plane direction
The spatial parameter, α = δz + jβ = −λ, can have a real part and an imaginary part. The
propagation constant is equal to the imaginary part, and the real part, δz, represents
the damping along the propagation direction.
• For a list of the studies available by physics interface, see The Wave
Optics Module Physics Interface Guide
• Studies and Solvers in the COMSOL Multiphysics Reference Manual
In this section:
DATA MANAGEMENT
With a very fine frequency step simulation, the solutions contain a lot of data. As a
result, the model file size will increase tremendously when it is saved. By selecting the
Store fields in output check box in the Values of Dependent Variables section of the
Frequency Domain study step settings, it is possible to define for what part of the model
the computed solution should be saved. When only S-parameters are of interest, it is
not necessary to store all of the field solutions. Instead, only store the field on the
selections for the port boundaries, as those will be used for the S-parameter
calculations.
In the Values of Dependent Variables section, change the selection in the Store fields in
output combo box from All to For selections and then add the explicit selections that
include the port boundaries. The explicit selection can be easily created from the port
feature by clicking Create Selection icon in the Boundary Selection settings once the
selection is specified.
DATA MANAGEMENT
The Store fields in output check box in the Values of Dependent Variables section can be
applied to the Frequency Domain, Modal study — if you are interested only in
S-parameters. By storing solutions only on port boundaries, the saved model file size
will decrease a lot.
ELECTROMAGNETIC QUANTITIES | 65
66 | CHAPTER 2: WAVE OPTICS MODELING
3
This chapter describes the physics interfaces found under the Optics>Wave Optics
branch ( ) when adding a physics interface.
In this chapter:
67
The Electromagnetic Waves,
Frequency Domain Interface
The Electromagnetic Waves, Frequency Domain (ewfd) interface ( ), found under the
Wave Optics branch ( ) when adding a physics interface, is used to solve for
time-harmonic electromagnetic field distributions.
For this physics interface, the maximum mesh element size should be limited to a
fraction of the wavelength. The domain size that can be simulated thus scales with the
amount of available computer memory and the wavelength. The physics interface
supports the Frequency Domain, Wavelength Domain, Eigenfrequency, Mode
Analysis, and Boundary Mode Analysis study types. The Frequency Domain and
Wavelength Domain study types are used for source driven simulations for a single
frequency/wavelength or a sequence of frequencies/wavelengths. The
Eigenfrequency study type is used to find resonance frequencies and their associated
eigenmodes in resonant cavities.
This physics interface solves the time-harmonic wave equation for the electric field.
When this physics interface is added, these default nodes are also added to the Model
Builder — Wave Equation, Electric, Perfect Electric Conductor, and Initial Values. Then,
from the Physics toolbar, add other nodes that implement, for example, boundary
conditions. You can also right-click Electromagnetic Waves, Frequency Domain to select
physics features from the context menu.
The Mode analysis study type is applicable only for 2D and 2D axisymmetric
cross-sections of waveguides and transmission lines where it is used to find allowed
propagating modes. Boundary mode analysis is used for the same purpose in 2D, 2D
axisymmetry, and 3D and applies to boundaries representing waveguide ports.
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
different physics interfaces, the name string must be unique. Only letters, numbers, and
underscores (_) are permitted in the Name field. The first character must be a letter.
SETTINGS
From the Formulation list, select whether to solve for the Full field (the default) or the
Scattered field.
For Scattered field select a Background wave type according to the following table:
TABLE 3-1: BACKGROUND WAVE TYPE BASED ON COMPONENT DIMENSION
User Defined
Enter the component expressions for the Background electric field Eb (SI unit: V/m).
The entered expressions must be differentiable.
Gaussian Beam
For Gaussian beam select the Gaussian beam type — Paraxial approximation (the default)
or Plane wave expansion.
A plane wave expansion with a finite number of plane waves included will make the
field periodic in the plane orthogonal to the main propagation direction. If the
separation between the transverse wave vector components, given by 2kt,max/(Nk-1),
is too small, replicas of the Gaussian beam background field can appear. To avoid that,
increase the value for the Wave vector count Nk.
The number of plane waves included in the expansion can be quite large, especially for
3D. For instance, using the default settings, 2*13*13 = 338 plane waves will be
included (the factor 2 accounts for the two possible polarizations for each wave
vector). Thus, initializing the plane-wave expansion for the Gaussian beam
background field can take some time in 3D.
For more information about the Gaussian beam theory, see Gaussian Beams as
Background Fields.
• Select a Beam orientation: Along the x-axis (the default), Along the y-axis, or for 3D
components, Along the z-axis.
• Enter a Beam radius w0 (SI unit: m). The default is 20π/ewfd.k0 m (10 vacuum
wavelengths).
• Enter a Focal plane along the axis p0 (SI unit: m). The default is 0 m.
• Enter the component expressions for the Background electric field amplitude,
Gaussian beam Ebg0 (SI unit: V/m).
• Enter a Wave number k (SI unit: rad/m). The default is ewfd.k0 rad/m. The wave
number must evaluate to a value that is the same for all the domains the scattered
field is applied to. Setting the Wave number k to a positive value, means that the wave
is propagating in the positive x-, y-, or z-axis direction, whereas setting the Wave
• Enter an Electric field amplitude E0 (SI unit: V/m). The default is 1 V/m.
• Enter a Roll angle (SI unit: rad), which is a right-handed rotation with respect to the
+x direction. The default is 0 rad, corresponding to polarization along the
+z direction.
• Enter a Pitch angle (SI unit: rad), which is a right-handed rotation with respect to
the +y direction. The default is 0 rad, corresponding to the initial direction of
propagation pointing in the +x direction.
Figure 3-1: Schematic of the directions for the wave vector k, the electric field E0, and the
roll, pitch and yaw rotations. The top image represents an initial wave propagating in the
x direction with a polarization along the z direction.
– jmϕ
E b ( r, ϕ, z ) = Ẽ b ( r, z )e
where
˜ ( r, z ) = E ( r̂ + jmϕ̂ )e –jkz ,
E b 0
and m is azimuthal mode number (+1 or -1) varying depending on the Circular
polarization type and Direction of propagation settings.
For 2D components, assign a wave vector component to the Out-of-plane wave number
field. For 2D axisymmetric components, assign an integer constant or an integer
parameter expression to the Azimuthal mode number field.
PHYSICS-CONTROLLED MESH
Select the Enable check box to use a physics-controlled mesh for the electromagnetic
problem. When selected, this invokes a parameter for the maximum mesh element size
in free space. The physics-controlled mesh automatically scales the maximum mesh
element size as the wavelength changes in different dielectric and magnetic regions. If
the model is configured by any periodic conditions, identical meshes are generated on
each pair of periodic boundaries. Perfectly matched layers are built with a structured
mesh, specifically, a swept mesh in 3D and a mapped mesh in 2D.
When Enable is selected, choose one of the four options for the Maximum mesh element
size control parameter — From study (the default), User defined, Frequency, or
Wavelength. When From study is selected, 1/5 of the vacuum wavelength from the
highest frequency defined in study step is used for the maximum mesh element size.
For the option User defined, enter a suitable Maximum element size in free space. For
example, 1/5 of the vacuum wavelength or smaller. When Frequency is selected, enter
the highest frequency intended to be used during the simulation. The maximum mesh
element size in free space is 1/5 of the vacuum wavelength for the entered frequency.
For the Wavelength option, enter the smallest vacuum wavelength intended to be used
during the simulation. The maximum mesh element size in free space is 1/5 of the
entered wavelength.
The maximum mesh element size in dielectric media is that in free space divided by the
square root of the product of the relative permittivity and permeability.
For Activate port sweep enter a Sweep parameter name to assign a specific name to the
parameter that controls the port number solved for during the sweep. Before making
the port sweep, the parameter must also have been added to the list of parameters in
the Parameters section of the Parameters node under the Global Definitions node.
For this physics interface, the S-parameters are subject to Touchstone file export. Click
Browse to locate the file, or enter a file name and path. Select an Output format:
Magnitude angle, Magnitude (dB) angle, or Real imaginary.
Enter a Reference impedance, Touchstone file export Zref (SI unit: Ω). The default is
50 Ω.
DEPENDENT VARIABLES
The dependent variables (field variables) are for the Electric field E and its components
(in the Electric field components fields). The name can be changed but the names of
fields and dependent variables must be unique within a model.
DISCRETIZATION
DOMAIN
BOUNDARY CONDITIONS
With no surface currents present the boundary conditions
n2 × ( E1 – E2 ) = 0
n2 × ( H1 – H2 ) = 0
need to be fulfilled. Because E is being solved for, the tangential component of the
electric field is always continuous, and thus the first condition is automatically fulfilled.
The second condition is equivalent to the natural boundary condition
–1 –1
– n × [ ( μ r ∇ × E ) 1 – ( μ r ∇ × E ) 2 ] = n × jωμ 0 ( H 1 – H2 ) = 0
In the COMSOL Multiphysics Reference Manual see Table 2-3 for links
to common sections and Table 2-4 to common feature nodes. You can
also search for information: press F1 to open the Help window or Ctrl+F1
to open the Documentation window.
2
∇ × ( μ r– 1 ∇ × E ) – k 0 ε rc E = 0
for the time-harmonic and eigenfrequency problems. The wave number of free space
k0 is defined as
ω
k 0 = ω ε 0 μ 0 = -----
c0
˜
E ( x, y, z ) = E ( x, y ) exp ( – ik z z ) .
˜ 2 ˜
( ∇ – ik z z ) × [ μ r–1 ( ∇ – ik z z ) × E ] – k 0 ε rc E = 0 ,
Notice that the ansatz above just explains how the wave equation is modified when the
out-of-plane wave vector component kz is not zero. As an example, for a plane wave
with a nonzero out-of-plane wave vector component, the electric field is of course
given by
˜
E ( x, y, z ) = E ( x, y ) exp ( – ik z z ) = A exp ( – i ( k x x + k y y + k z z ) ) ,
where A is a constant amplitude and kx, ky, and kz are the wave vector components.
In 2D axisymmetry, the electric field varies with the azimuthal mode number m as
˜
E ( r, ϕ, z ) = E ( r, z ) exp ( – imϕ ) .
∇ – im ˜ ˜
----- ϕ × μ r–1 ∇ – i m
----- ϕ × E – k 0 ε rc E = 0 ,
2
r r
ω
Q fact = ---------
2δ
Using the relation εr = n2, where n is the refractive index, the equation can
alternatively be written
2 2
∇ × ( ∇ × E ) – k0 n E = 0
When the equation is written using the refractive index, the assumption is that μr = 1
and σ = 0 and only the constitutive relations for linear materials are available. When
solving for the scattered field the same equations are used but E = Esc + Ei and Esc is
the dependent variable.
Relative Permittivity
When Relative permittivity is selected, the default Relative permittivity εr takes values
From material. For User defined select Isotropic, Diagonal, Symmetric, or Anisotropic and
enter values or expressions in the field or matrix.
Refractive Index
When Refractive index is selected, the default Refractive index n and Refractive index,
imaginary part k take the values From material. To specify the real and imaginary parts
of the refractive index and assume a relative permeability of unity and zero
Dielectric Loss
When Dielectric loss is selected, the default Relative permittivity ε′ and Relative
permittivity (imaginary part) ε″ take values From material. For User defined select
M 2
fj ωP
εr ( ω ) = ε∞ + ---------------------------------------
2
ω 0j – ω + iΓ j ω
2
-
j=1
When Drude-Lorentz dispersion model is selected, the default Relative permittivity, high
frequency ε∞ (dimensionless) takes its value From material. For User defined select
Isotropic, Diagonal, Symmetric, or Anisotropic and enter a value or expression in the field
or matrix.
In the table, enter values or expressions in the columns for the Oscillator strength,
Resonance frequency (rad/s), and Damping in time (rad/s).
Δε k
ε ( ω ) = ε∞ + ---------------------
1 + iωτ k
-
k
When Debye dispersion model is selected, the default Relative permittivity, high
frequency ε∞ (dimensionless) takes its value From material. For User defined select
Isotropic, Diagonal, Symmetric, or Anisotropic and enter a value or expression in the field
or matrix.
2
Bk λ
2
n (λ) = 1 + ------------------
-
2
k λ – Ck
When Sellmeier dispersion model is selected, in the table, enter values or expressions in
the columns for B and C (m^2).
MAGNETIC FIELD
Select the Constitutive relation — Relative permeability (the default) or Magnetic losses.
• For Relative permeability the relative permeability μr uses values From material. For
User defined select Isotropic, Diagonal, Symmetric, or Anisotropic based on the
characteristics of the magnetic field, and then enter values or expressions in the field
or matrix.
• For Magnetic losses the default values for Relative permeability (real part) μ′ and
Relative permeability (imaginary part) μ″ are taken From material. For User defined
enter different values.
CONDUCTION CURRENT
By default, the Electrical conductivity σ (SI unit: S/m) uses values From material.
• For User defined select Isotropic, Diagonal, Symmetric, or Anisotropic based on the
characteristics of the current and enter values or expressions in the field or matrix.
• For Linearized resistivity the default values for the Reference temperature Tref (SI
unit: K), Resistivity temperature coefficient α (SI unit: 1/K), and Reference resistivity
ρ0 (SI unit: Ω⋅m) are taken From material. For User defined enter other values or
expressions for any of these variables.
Initial Values
The Initial Values node adds an initial value for the electric field that can serve as an
initial guess for a nonlinear solver. Add additional Initial Values nodes from the Physics
toolbar.
INITIAL VALUES
Enter values or expressions for the initial values of the components of the Electric field
E (SI unit: V/m). The default values are 0 V/m.
J = σE + J e
Far-Field Domain
To set up a far-field calculation, add a Far-Field Domain node and specify the far-field
domains in its Settings window. Use Far-Field Calculation subnodes (one is added by
default) to specify all other settings needed to define the far-field calculation. If a
Perfectly Matched Layer (PML) node has been added before adding the Far-Field
Domain, all of the domains in the Electromagnetic Waves, Frequency Domain interface
adjacent to the PML are automatically selected by default. If there is no PML, all of
the domains are selected. The selection can be modified. In that case, select only a
Far-Field Calculation
A Far-Field Calculation subnode is added by default to the Far-Field Domain node and
is used to select boundaries corresponding to a single closed surface surrounding all
radiating and scattering objects. By default, all exterior boundaries of the Far-Field
Domain are selected. If a Perfectly Matched Layer (PML) node has been added before
adding the Far-Field Domain, all exterior boundaries of the Far-Field Domain adjacent
to the PML are selected. Symmetry reduction of the geometry makes it relevant to
select boundaries defining a nonclosed surface. Also use this feature to indicate
symmetry planes and symmetry cuts applied to the geometry, and whether the selected
boundaries are defining the inside or outside of the far field domain; that is, to say
whether they are facing away from infinity or toward infinity.
FAR-FIELD CALCULATION
Enter a Far-field variable name. The default is Efar.
Select as needed the Symmetry in the x=0 plane, Symmetry in the y=0 plane, or
Symmetry in the z=0 plane check boxes to use it your model when calculating the
far-field variable. The symmetry planes have to coincide with one of the Cartesian
coordinate planes.
When a check box is selected, also choose the type of symmetry to use from the
Symmetry type list that appears — Symmetry in E (PMC) or Symmetry in H (PEC). The
selection should match the boundary condition used for the symmetry boundary.
Using these settings, include the parts of the geometry that are not in the model for
symmetry reasons in the far-field analysis.
If perfectly matched layers are added to the model after the Far-Field
Domain is configured, then it is necessary to press the Reset Far-Field
Boundaries button to reassign all exterior boundaries.
Polarization
The Polarization node adds an externally generated polarization Pi, which contributes
to the total polarization
P = ε 0 ( ε r – 1 )E + Pi .
i
Add Polarization nodes from the Physics toolbar or by right-clicking the physics
interface and selecting the Polarization item from the context menu.
POLARIZATION
Enter the components (x, y, and z for 3D components for example) of the Polarization
Pi (SI unit: C/m2).
is a special case of the electric field boundary condition that sets the tangential
component of the electric field to zero. It is used for modeling of a lossless metallic
surface (for example, a ground plane) or as a symmetry type boundary condition. It
imposes symmetry for magnetic fields and “magnetic currents” and antisymmetry for
electric fields and electric currents. It supports induced electric surface currents and
thus any prescribed or induced electric currents (volume, surface, or edge currents)
flowing into a perfect electric conductor boundary is automatically balanced by
induced surface currents.
Js
J
I'
I
Js
Figure 3-2: The perfect electric conductor boundary condition is used on exterior and
interior boundaries representing the surface of a lossless metallic conductor or (on exterior
boundaries) representing a symmetry cut. The shaded (metallic) region is not part of the
model but still carries effective mirror images of the sources. Note also that any current
flowing into the boundary is perfectly balanced by induced surface currents. The
tangential electric field vanishes at the boundary.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
For information about the Constraint Settings section, see Constraint Settings in the
COMSOL Multiphysics Reference Manual.
n×H = 0
Js=0
I'
I
J=0
Figure 3-3: The perfect magnetic conductor boundary condition is used on exterior
boundaries representing the surface of a high impedance region or a symmetry cut. The
shaded (high impedance) region is not part of the model but nevertheless carries effective
mirror images of the sources. Note also that any electric current flowing into the boundary
is forbidden as it cannot be balanced by induced electric surface currents. The tangential
magnetic field vanishes at the boundary. On interior boundaries, the perfect magnetic
conductor boundary condition literally sets the tangential magnetic field to zero which in
addition to setting the surface current density to zero also makes the tangential electric
field (and in dynamics the tangential electric field) discontinuous.
Port
Use the Port node where electromagnetic energy enters or exits the model. A port can
launch and absorb specific modes. Use the boundary condition to specify wave type
In 3D, also right-click these subnodes are available from the context menu (right-click
the parent node) or from the Physics toolbar, Attributes menu:
• Circular Port Reference Axis to determine a reference direction for the modes. This
subnode is selected from the Points submenu when Circular is selected as the type of
port.
• Periodic Port Reference Point to uniquely determine reciprocal lattice vectors. This
subnode is selected from the Points submenu when Periodic is selected as the type of
port.
PORT PROPERTIES
Enter a unique Port name. Only nonnegative integer numbers can be used as Port name
as it is used to define the elements of the S-parameter matrix and numeric port names
are also required for port sweeps and Touchstone file export.
Select the Type of Port — User defined, Numeric, Rectangular, Coaxial, Circular, or
Periodic.
Periodic ports are available in 3D and 2D. Circular and Coaxial ports are available in 3D
and 2D axisymmetry.
Numeric ports require a Boundary Mode Analysis study type. It should appear before the
frequency or wavelength domain study node in the study branch of the model tree. If
more than one numeric port is needed, use one Boundary Mode Analysis node per
port and assign each to the appropriate port. Then, it is best to add all the studies;
Boundary Mode Analysis 1, Boundary Mode Analysis 2,..., Frequency Domain 1 (or
Wavelength Domain 1), manually. Numeric ports are by default computed for the
deformed mesh whereas other types of ports compute the mode shape using geometry
information.
Then select a Slit type — PEC-backed (the default) or Domain-backed. The PEC-backed
type makes the port on interior boundaries perform as it does on exterior boundaries.
The Domain-backed type can be combined with perfectly matched layers to absorb the
excited mode from a source port and other higher order modes.
Click Toggle Power Flow Direction button to define the power flow for the port. For an
excited port, the power flow should point in to the excited domain and for a listener
port the power flow should point out from the excited domain. The power flow
direction is visualized with a red arrow on the port boundary in the Graphics window.
User Defined
For User defined specify the eigenmode of the incoming wave at the port. Even if Wave
excitation at this port is set to Off, the mode field should represent the incoming wave
that corresponds to the actual outgoing wave. The mode field can be entered with an
arbitrary amplitude and is normalized internally.
• Enter the components of the Electric mode field E0 (SI unit: V/m) or the Magnetic
mode field H0 (SI unit: A/m). The entered expressions must be differentiable.
• Enter the Propagation constant β (SI unit: rad/m). This is frequency dependent for
all but TEM modes and a correct frequency-dependent expression must be used.
Rectangular
For Rectangular specify a unique rectangular mode.
For 2D components, to excite the fundamental mode, select the mode type Transverse
electromagnetic (TEM), since the rectangular port represents a parallel-plate waveguide
port that can support a TEM mode. Only TE modes are possible when solving for the
out-of-plane vector component, and only TM and TEM modes are possible when
solving for the in-plane vector components. There is only a single mode number,
which is selected from a list.
Coaxial
In 2D axisymmetry, Coaxial does not support nonzero azimuthal mode number. The
Azimuthal mode number in the Physics interface should be defined as zero.
For 3D components, enter the Mode number, for example, 11 for a TE11 mode, or 01
for a TM01 mode. When Circular is selected as the type of port in 3D, the Circular Port
Reference Axis subnode is available from the context menu (right-click the parent
node) or from the Physics toolbar, Attributes menu. It defines the orientation of fields
on a port boundary.
Periodic
For Periodic specify parameters for the mode field. When Periodic is selected, the
Diffraction Order port subnode is available from the context menu (right-click the
parent node) or from the Physics toolbar, Attributes menu.
Select a Input quantity — Electric field or Magnetic field and define the mode field
amplitude for the incoming wave at the port. Even if Wave excitation at this port is set
to Off, the mode field amplitude should represent the incoming wave that corresponds
to the actual outgoing wave.
• For 2D components and if the Input quantity is set to Electric field, define the Electric
mode field amplitude. For example, for a TE wave set the x, y, and z components to
0, 0, 1. Similarly, if the Input quantity is set to Magnetic field, define the Magnetic
mode field amplitude. For a TM wave set the x, y, and z components to 0, 0, and 1.
• Define the Angle of incidence, if Wave excitation at this port is On.
For 3D components, if Wave excitation at this port is On, define the Elevation angle of
incidence and Azimuth angle of incidence. The Elevation angle of incidence α1 and Azimuth
angle of incidence α2 are used in the relations
k = k parellel + k perpendicular
The Azimuth angle of incidence is the counterclockwise rotating angle from the
primitive vector a1 around the axis built with Periodic Port Reference Point and n.
For periodic ports with hexagonal port boundaries, the definition of the vector a1 is
slightly different from the default definition. In this case, the unit cell is actually a
rhomboid, with primitive vectors pointing in other directions than the side vectors of
the hexagon. Thus, for a hexagonal periodic port, the vector a1 is defined along one
of the sides of the hexagon, and it is not one of the primitive vectors of the hexagonal
point lattice. The Azimuth angle of incidence α2 is still measured from the vector a1,
even though this vector now refers to a side vector of the hexagonal port boundary and
not a primitive vector.
For 2D components define the Angle of incidence. The Angle of incidence α is defined
by the relation
k × n = k sin αz
where k is the projection of the wave vector in the xy-plane, n is the normalized
normal vector to the boundary, k is the magnitude of the projected wave vector in the
xy-plane, and z is the unit vector in the z direction.
Notice that the mode field defined for the Periodic port assumes
homogeneous isotropic material properties in the domain adjacent to the
selected port boundary.
• Select the Include in automatic diffraction order calculation check box to add
Diffraction Order subnodes to the selected Periodic port, when the Compute
Diffraction Order button is clicked from the exciting Periodic port.
• Define the Refractive index, real part at the boundary.
• Define the Maximum frequency — From study (the default) or User defined. When
From study is selected, the Maximum frequency is taken from the study step associated
with the physics interface. For User defined, enter the maximum frequency fmax (SI
unit: Hz). The default value is 0 Hz. If a single frequency is used, insert the
frequency, or if a frequency sweep is performed, insert the maximum frequency of
the sweep. This parameter is only available when Wave excitation at this port is On.
When all parameters are defined, click the Compute Diffraction Orders button from the
exciting Periodic port to automatically create Diffraction Order ports as subnodes to all
Periodic ports having the Include in automatic diffraction order calculation check box
selected.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
For information about the Constraint Settings section, see Constraint Settings in the
COMSOL Multiphysics Reference Manual.
Diffraction Order
The Diffraction Order port is available in 3D and 2D. When the Type of Port is set to
Periodic under Port Properties, this subnode is available from the context menu
(right-click the Port parent node) or from the Physics toolbar, Attributes menu.
Use the Diffraction Order port to define diffraction orders from a periodic structure.
Normally a Diffraction Order node is added automatically during the Periodic port
setup. Additional Diffraction Order ports subnodes are available from the context menu
(right-click the parent node) or from the Physics toolbar, Attributes menu.
PORT PROPERTIES
Enter a unique Port name. Only nonnegative integer numbers can be used as Port name
as it is used to define the elements of the S-parameter matrix and numeric port names
are also required for port sweeps and Touchstone file export.
Diffraction Order
Specify an integer constant or an integer parameter expression for the Diffraction order
m (the default is 0) and in 3D n (the default is 0).
Note that In-plane vector and Out-of-plane vector are based on the plane of diffraction
which is constructed with the diffraction wave vector and the outward normal vector
of the port boundary. The diffraction wave vector is defined by
k diffraction,parallel = k F + MG 1 + NG 2
2 2
k diffraction,perpendicular = k – k diffraction,parallel
In-plane vector lies on the plane of diffraction while Out-of-plane vector is normal to the
plane of diffraction.
For a 2D component, In-plane vector is available when the settings for the physics
interface is set to either In-plane vector or Three-component vector under Electric Field
Components Solved For. Out-of-plane vector is available when the settings for the
physics interface is set to either Out-of-plane vector or Three-component vector under
Electric Field Components Solved For.
Enter a value or expression for the Mode phase θin (SI unit: rad). The default is
0 radians. The Mode phase setting is further discussed for the Port feature.
Notice that the mode field defined for the Diffraction Order feature
assumes isotropic material properties in the domain adjacent to the
selected feature boundary.
The Periodic Port Reference Point is used to uniquely identify two primitive unit cell
vectors, a1 and a2, and two reciprocal lattice vectors, G1 and G2. These reciprocal
vectors are defined in terms of the unit cell vectors, a1 and a2, tangent to the edges
shared between the port and the adjacent periodic boundary conditions. G1 and G2
are defined by the relation
a1 × a2
---------------------
- = n
a1 × a2
a2 × n n × a1
G 1 = 2π --------------------------- and G 2 = 2π ---------------------------
a1 ⋅ a2 × n a1 ⋅ a2 × n
where n is the outward unit normal vector to the port boundary. If there are multiple
points defined in the selection list, only the last point is used.
POINT SELECTION
The primitive unit cell vectors, a1 and a2 are defined from two edges sharing the
Periodic Port Reference Point on a port boundary. The two vectors can have unequal
lengths and are not necessarily orthogonal. They start from the Periodic Port Reference
Point.
For listener (passive, observation, and not excited) ports, if the outward normal vector
on the listener port boundary is opposite to that of the source port, the listener port
reference point needs to be mirrored from the source port reference point based on
the center coordinate of the model domain. For example, if the source port reference
point is at {-1,-1,1} in a cubic domain around the origin, the mirrored listener port
If the lattice vectors are collinear with two Cartesian axes, then the lattice vectors can
be defined without the Periodic Port Reference Point. For the port where n points along
a positive Cartesian direction, a1 and a2 are also assigned to point along positive
Cartesian directions. Conversely, for the port where n points along a negative
Cartesian direction, a1 and a2 are assigned to point along negative Cartesian
directions. The condition a1 × a2 || n is true on both ports. For example, if n = z, then
a1/|a1| = x and a2/|a2| = y and if n = −z, then a1/|a1| = −x and a2/|a2| = −y.
Electric Field
The Electric Field boundary condition
n × E = n × E0
specifies the tangential component of the electric field. It should in general not be used
to excite a model. Consider using the Port or Scattering Boundary Condition instead.
It is provided mainly for completeness and for advanced users who can recognize the
special modeling situations when it is appropriate to use. The commonly used special
case of zero tangential electric field is described in the Perfect Electric Conductor
section.
ELECTRIC FIELD
Enter the value or expression for the components of the Electric field E0
(SI unit: V/m).
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
For information about the Constraint Settings section, see Constraint Settings in the
COMSOL Multiphysics Reference Manual.
n × H = n × H0
MAGNETIC FIELD
Enter the value or expression for the components of the Magnetic field H0
(SI unit: A/m).
– jk ( n ⋅ r ) – jk ( k ⋅ r )
E = E sc e + E0 e Plane scattered wave
– jk ( n ⋅ r )
e – jk ( k ⋅ r )
E = E sc ------------------------ + E 0 e Cylindrical scattered wave
r
– jk ( n ⋅ r )
e – jk ( k ⋅ r )
E = E sc ------------------------ + E 0 e Spherical scattered wave
rs
The field E0 is the incident plane wave that travels in the direction k. The boundary
condition is transparent for incoming (but not outgoing) plane waves with any angle
of incidence.
– jk ( k ⋅ ( r – r ref ) )
E0 e ,
where rref is a reference point determined as the average point from the point selection
in the Reference Point subnode.
The boundary is only perfectly transparent for scattered (outgoing) waves of the
selected type at normal incidence to the boundary. That is, a plane wave at oblique
• For cylindrical waves, specify around which cylinder axis the waves are cylindrical.
Do this by specifying one point at the cylinder axis and the axis direction.
• For spherical waves, specify the center of the sphere around which the wave is
spherical.
The domain material adjacent to the boundary where the Scattering Boundary
Condition is applied can be lossy.
If the problem is solved for the eigenfrequency or the scattered field, the boundary
condition does not include the incident wave.
– jk ( n ⋅ r )
E sc = E sc e Plane scattered wave
– jk ( n ⋅ r )
e
E sc = E sc ------------------------ Cylindrical scattered wave
r
– jk ( n ⋅ r )
e
E sc = E sc ------------------------ Spherical scattered wave
rs
If the Incident field is not set to No incident field, edit the Incident wave direction kdir
for the vector coordinates. The default direction is in the opposite direction to the
boundary normal. For 2D axisymmetry, the Incident wave direction kdir should be
parallel or anti-parallel to the symmetry axis.
Select a Scattered wave type for which the boundary is absorbing — Plane wave (the
default), Spherical wave, or Cylindrical wave.
• For the Electromagnetic Waves, Frequency Domain interface, select an Order —First
order (the default) or Second order.
• For Cylindrical wave also enter coordinates for the Source point r0 (SI unit: m) and
Source axis direction raxis (dimensionless). For 2D the Source axis direction is
Select the Dispersion and absorption model that will be used when calculating the wave
number and attenuation constant for the incident and scattered waves — Low loss
When the Dispersion and absorption model is set to Low loss approximation
the refractive index is calculated from the relative permittivity and the
relative permeability as
n = εr μr .
1 μ0 μr 1
γ = --- σ -----------
- = --- σZ c ,
2 ε0 εr 2
When the Dispersion and absorption model is set to High loss, the real and
the imaginary parts of the complex refractive index is solved for from the
real and the imaginary parts of the relative permittivity, using the relations
2 2
n – κ = ε' r μ r
and
σμ r
2nκ = ε'' r μ r = --------- .
ωε 0
ω
γ = ---- κ .
c
Reference Point
The Reference Point subnode is available only when there is an available incident field
defined in the parent node. Then this subnode is available from the context menu
(right-click the Port parent node) or from the Physics toolbar, Attributes menu.
The Reference Point subnode defines a reference position rref that is calculated as the
average position from the point selection in the Reference Point subnode.
In the parent node, the incident field is then defined using the reference position:
POINT SELECTION
Select the points that should be used when calculating the reference position. The
reference position is calculated as the average position of the selected points.
μ0 μr
------------ n × H + E – ( n ⋅ E )n = ( n ⋅ E s )n – E s
εc
is used at boundaries where the field is known to penetrate only a short distance
outside the boundary. This penetration is approximated by a boundary condition to
avoid the need to include another domain in the model. Although the equation is
identical to the one in the low-reflecting boundary condition, it has a different
interpretation. The material properties are for the domain outside the boundary and
not inside, as for low-reflecting boundaries. A requirement for this boundary condition
to be a valid approximation is that the magnitude of the complex refractive index
με c
N = -----------
-
μ1 ε1
where μ1 and ε1 are the material properties of the inner domain, is large; that is,
|N| >> 1.
When used with the Electromagnetic Waves, Beam Envelopes interface, the
propagation direction in the exterior layer can be specified. Setting the propagation
direction to be in the normal direction is the default option and results in the behavior
as described above. However, setting the propagation direction to be given by the
wave vector direction takes the tangential wave vector components from the Wave
Vectors settings for the physics and the longitudinal component is derived to make the
wave number satisfy the wave number for the exterior layer. This option is useful if the
exterior domain is a dielectric material.
The source electric field Es can be used to specify a source surface current on the
boundary.
Js
PROPAGATION DIRECTION
This section is only available for the Electromagnetic Waves, Beam Envelopes
interface. Select a Propagation direction — Normal direction (the default) or From wave
vector. The Normal direction option assumes that the wave in the exterior material
propagates essentially in the normal direction, whereas the From wave vector option
assumes that the tangential wave vector component is continuous at the boundary, as
specified by the wave vectors k1 and k2 for the Electromagnetic Waves, Beam
Envelopes interface. The normal component for the wave vector in the exterior
material is obtained from the wave number, given the material parameters of the
exterior domain. Thus, this option implements Snell’s law of refraction at the
boundary, which makes this option useful also for dielectric exterior materials.
–n × H = Js
n2 × ( H1 – H2 ) = Js
specifies a surface current density at both exterior and interior boundaries, respectively.
The current density is specified as a three-dimensional vector, but because it needs to
flow along the boundary surface, COMSOL Multiphysics projects it onto the
boundary surface and neglects its normal component. This makes it easier to specify
the current density and avoids unexpected results when a current density with a
component normal to the surface is given.
For the Surface Current Density subnode, select Side — Upside (the default) or Downside
to define on which side the Surface Current Density is applied. The red arrow visualized
on the selected boundaries always indicates the upside.
n × E = J ms
n 2 × ( E 1 – E 2 ) = – J ms
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
For information about the Constraint Settings section, see Constraint Settings in the
COMSOL Multiphysics Reference Manual.
( Z S E t1 – Z T E t2 )
J s1 = --------------------------------------------
2 2
-
ZS – ZT
( Z S E t2 – Z T E t1 )
J s2 = --------------------------------------------
2 2
-
ZS – ZT
– jωμ 1
Z S = ------------- ----------------------
k tan ( kd )
– jωμ 1
Z T = ------------- ---------------------
k sin ( kd )
k = ω ( ε + ( σ ⁄ ( jω ) ) )μ
The thickness of the layer should also be less than the radius of curvature
for the boundary.
PROPAGATION DIRECTION
This section is only available for the Electromagnetic Waves, Beam Envelopes
interface. Select a Propagation direction — Normal direction (the default) or From wave
vector. The Normal direction option assumes that the waves in the layer propagate
essentially in the normal direction, whereas the From wave vector option assumes that
the tangential wave vector component is continuous at the layer boundaries, as
specified by the wave vectors k1 and k2 for the Electromagnetic Waves, Beam
Envelopes interface. The normal component for the wave vector in the layer is
obtained from the wave number, given the specified material parameters. Thus, this
option implements Snell’s law of refraction for the layer, which makes this option
useful also for dielectric layers.
Select the Electrically thick layer check box (unselected by default) to make the two
domains adjacent to the boundary uncoupled. When the Electrically thick layer check
box is unselected, enter a Thickness d (SI unit: m). The default is 0.01 m.
Periodic Condition
The Periodic Condition sets up a periodicity between the selected boundaries. The
Destination Selection subnode is available from the context menu (right-click the parent
node) or from the Physics toolbar, Attributes menu.
BOUNDARY SELECTION
The software automatically identifies the boundaries as either source boundaries or
destination boundaries This works fine for cases like opposing parallel boundaries. To
control the destination, add a Destination Selection subnode. By default it contains the
selection that COMSOL Multiphysics has identified.
• Continuity to make the electric field periodic (equal on the source and destination),
• Antiperiodicity to make it antiperiodic, or
• Floquet periodicity (The Electromagnetic Waves, Frequency Domain Interface and
The Electromagnetic Waves, Beam Envelopes Interface only) to use a Floquet
periodicity (Bloch-Floquet periodicity).
- For Floquet periodicity also enter the source for the k-vector for Floquet periodicity.
- For User defined specify the components of the k-vector for Floquet periodicity kF
(SI unit: rad/m).
- For From periodic port the k-vector for Floquet periodicity kF is obtained from the
Periodic Port settings.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
For information about the Constraint Settings section, see Constraint Settings in the
COMSOL Multiphysics Reference Manual.
ORIENTATION OF SOURCE
For information about the Orientation of Source section, see Orientation of Source and
Destination in the COMSOL Multiphysics Reference Manual.
Magnetic Current
The Magnetic Current node specifies a magnetic line current along one or more edges.
For a single Magnetic Current source, the electric field is orthogonal to both the line
MAGNETIC CURRENT
Enter a value for the Magnetic current Im (SI unit: V).
Edge Current
The Edge Current node specifies an electric line current along one or more edges.
EDGE CURRENT
Enter an Edge current I0 (SI unit: A).
DIPOLE SPECIFICATION
Select a Dipole specification — Magnitude and direction or Dipole moment.
DIPOLE PARAMETERS
Based on the Dipole specification selection:
• For Magnitude and direction enter coordinates for the Electric current dipole moment
direction np and Electric current dipole moment, magnitude p (SI unit: A·m).
• For Dipole moment enter coordinates for the Electric current dipole moment p (SI
unit: A·m).
DIPOLE SPECIFICATION
Select a Dipole specification — Magnitude and direction or Dipole moment.
• For Magnitude and direction enter coordinates for the Magnetic dipole moment
direction nm and Magnetic dipole moment, magnitude m (SI unit: m2·A).
• For Dipole moment enter coordinates for the Magnetic dipole moment m (SI unit:
m2·A).
When this physics interface is added, these default nodes are also added to the Model
Builder — Wave Equation, Electric, Perfect Electric Conductor, and Initial Values. Then,
from the Physics toolbar, add other nodes that implement, for example, boundary
conditions and mass sources. You can also right-click Electromagnetic Waves, Transient
to select physics features from the context menu.
Except where indicated, most of the settings are the same as for The Electromagnetic
Waves, Frequency Domain Interface.
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
different physics interfaces, the name string must be unique. Only letters, numbers, and
underscores (_) are permitted in the Name field. The first character must be a letter.
The default Name (for the first physics interface in the model) is ewt.
COMPONENTS
This section is available for 2D and 2D axisymmetric components.
DEPENDENT VARIABLES
The dependent variable (field variable) is for the Magnetic vector potential A. The name
can be changed but the names of fields and dependent variables must be unique within
a model.
DISCRETIZATION
The domain, boundary, edge, point, and pair nodes are available from the Physics
ribbon toolbar (Windows users), Physics context menu (Mac or Linux users), or
right-click to access the context menu (all users).
DOMAIN
These nodes are unique for this physics interface and described in this section:
BOUNDARY CONDITIONS
With no surface currents present the boundary conditions
n2 × ( E1 – E2 ) = 0
n2 × ( H1 – H2 ) = 0
need to be fulfilled. Depending on the field being solved for, it is necessary to analyze
these conditions differently. When solving for A, the first condition can be formulated
in the following way.
∂A 2 ∂A 1 ∂
n2 × ( E1 – E2 ) = n2 × – = ( n2 × ( A2 – A1 ) )
∂t ∂t ∂t
The tangential component of the magnetic vector potential is always continuous and
thus the first condition is fulfilled. The second condition is equivalent to the natural
boundary condition.
–1 –1 –1
–n × ( μr ∇ × A1 – μr ∇ × A2 ) = –n × μr ( H1 – H2 ) = 0
These nodes and subnodes are available and described for the Electromagnetic Waves,
Frequency Domain interface (listed in alphabetical order):
In the COMSOL Multiphysics Reference Manual see Table 2-3 for links
to common sections and Table 2-4 to common feature nodes. You can
also search for information: press F1 to open the Help window or Ctrl+F1
to open the Documentation window.
∂A ∂ ∂A
+ μ0 ε0 εr + ∇ × ( μr ∇ × A ) = 0
–1
μ0 σ
∂t ∂t ∂ t
for transient problems with the constitutive relations B = μ0μrH and D = ε0εrE. Other
constitutive relations can also be handled for transient problems.
Refractive Index
When Refractive index is selected, the default Refractive index n (dimensionless) takes
the value From material. To specify the refractive index and assume a relative
permeability of unity and zero conductivity, for one or both of the options, select User
defined then choose Isotropic, Diagonal, Symmetric, or Anisotropic. Enter values or
expressions in the field or matrix.
Notice that only the real part of the refractive index is used for the
transient formulation.
Polarization
For Polarization enter coordinates for the Polarization P (SI unit: C/m2).
D = ε0 ε∞ E + Pn ,
n=1
MAGNETIC FIELD
This section is available if Relative permittivity, Polarization, or Remanent electric
displacement are chosen as the Electric displacement field model.
Select the Constitutive relation — Relative permeability (the default), Remanent flux
density, or Magnetization.
Relative Permeability
For Relative permeability the relative permeability μr uses values From material. For User
defined select Isotropic, Diagonal, Symmetric, or Anisotropic based on the characteristics
of the magnetic field, and then enter values or expressions in the field or matrix.
Magnetization
For Magnetization enter coordinates for M (SI unit: A/m).
CONDUCTION CURRENT
This section is available if Relative permittivity, Polarization, or Remanent electric
displacement are chosen as the Electric displacement field model.
By default, the Electrical conductivity σ (SI unit: S/m) uses values From material.
• For User defined select Isotropic, Diagonal, Symmetric, or Anisotropic based on the
characteristics of the current and enter values or expressions in the field or matrix.
• For Linearized resistivity the default values for the Reference temperature Tref (SI
unit: K), Resistivity temperature coefficient α (SI unit: 1/K), and Reference resistivity
ρ0 (SI unit: Ωm) use values From material. For User defined enter other values or
expressions for any of these variables.
INITIAL VALUES
Enter values or expressions for the initial values of the components of the magnetic
vector potential A (SI unit: Wb/m) and its time derivative ∂A/∂t (SI unit: V/m). The
default values are 0 Wb/m and 0 V/m, respectively.
Drude-Lorentz Polarization
This subfeature is available only when Drude-Lorentz Dispersion Model is selected as the
Electric displacement field model in the Wave Equation, Electric feature node. Then the
subnodes are made available from the context menu (right-click the parent node) as
well as from the Physics toolbar, Attributes menu.
D = ε0 ε∞ E + Pn ,
n=1
∂2 ∂ 2 2
-------2- + Γ n ----- + ω n P n = ε 0 f n ω p E .
∂t ∂t
Enter values or expressions for the Oscillator strength fn (SI unit: 1), the Resonance
frequency ωn (SI unit: rad/s), and the Damping in time coefficient Γn (SI unit: rad/s).
INITIAL VALUES
Enter values or expressions for the initial values of the components of the
Drude-Lorentz polarization Pn (SI unit: C/m2) and its time derivative ∂Pn/∂t (SI
unit: A/m2). The default values are 0 C/m2 and 0 A/m2, respectively.
This physics interface solves two first-order partial differential equations (Faraday’s law
and Maxwell-Ampère’s law) for the electric and magnetic fields using the time explicit
discontinuous Galerkin method.
When this physics interface is added, these default nodes are also added to the Model
Builder — Wave Equations, Perfect Electric Conductor, and Initial Values. Then, from the
Physics toolbar, add other nodes that implement, for example, boundary conditions.
You can also right-click Electromagnetic Waves, Time Explicit to select physics features
from the context menu.
The interface includes absorbing layers that are used to set up effective nonreflecting
like boundary conditions. These features are added from the Definitions toolbar, by
clicking Absorbing Layer. If COMSOL Multiphysics is not running in full-screen mode
nor in a large window, Absorbing Layer is accessible on the Definitions toolbar by first
clicking Coordinate Systems and then Absorbing Layer. You can also right-click
Definitions in the Model Builder and select Absorbing Layer from the context menu.
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
different physics interfaces, the name string must be unique. Only letters, numbers, and
underscores (_) are permitted in the Name field. The first character must be a letter.
The default Name (for the first physics interface in the model) is teew.
COMPONENTS
This section is available for 2D and 2D axisymmetric components.
• Full wave (the default) to solve using a full three-component vector for the electric
field E and the magnetic field H.
• E in plane (TM wave) to solve for the electric field vector components in the modeling
plane and one magnetic field vector component perpendicular to the plane,
assuming that there is no electric field perpendicular to the plane and no magnetic
field components in the plane.
• H in plane (TE wave) to solve for the magnetic field vector components in the
modeling plane and one electric field vector component perpendicular to the plane.
DEPENDENT VARIABLES
The dependent variables (field variables) are for the Electric field vector E and for the
Magnetic field vector H. The name can be changed but the names of fields and
dependent variables must be unique within a model.
DISCRETIZATION
Wave Equations
The Wave Equations node is the main node for the Electromagnetic Waves, Time
Explicit interface. The governing transient equations can be written in the form
∂D
∇ × H = σE + -------
∂t
∂B
∇ × E = – -------
∂t
∂E
ε 0 ε r ------- – ∇ × H + σE = 0
∂t
∂H
μ 0 μ r -------- + ∇ × E = 0
∂t
MATERIAL PROPERTIES
The default Relative permittivity εr (dimensionless), Relative permeability μr
(dimensionless), and Electrical conductivity σ (SI unit: S/m) take values From material.
NUMERICAL PARAMETERS
The defaults for each parameter are as follows:
• Lax-Friedrichs flux parameter for E field τE (SI unit: S), the default is 0.5/Z for
Ampere’s law.
• Lax-Friedrichs flux parameter for H fieldτH (SI unit:Ω), the default is 0.5 Z for
Faraday’s law, where Z is the impedance of vacuum.
• Estimate of maximum wave speed cmax (SI unit: m/s) the default is taken from the
speed of light in a vacuum c_const.
FILTER PARAMETERS
The filter provides higher-order smoothing of nodal discontinuous Galerkin
formulations and is intended to be used for absorbing layers, but you can also use it to
stabilize linear wave problems with highly varying coefficients. The filter is constructed
by transforming the solution (in each global time step) to an orthogonal polynomial
representation, multiplying with a damping factor and then transforming back to the
(Lagrange) nodal basis. Select the Activate check box to use this filter.
–1
VΛV
1, 0 ≤ η ≤ η c
Λ mm = σ(η) = η – η c 2s
--------------
-
– α 1 – η c
e , ηc ≤ η ≤ 1
where
im
η = η ( m ) = -------
Np
and Np is the basis function and im the polynomial order for coefficient m. α (default
value: 36), ηc (default value: 1), and s (default value: 3) are the filter parameters that
When using Absorbing Layer features, add an additional Wave Equations feature for the
corresponding domain selection. Select the Activate check box and add filter
parameters. An example of a filter parameter combination that can be used for a Wave
Equations feature active on an Absorbing Layer domain selection is α = 0.1, ηc = 0.01,
and s = 2. However, other combinations could work better, depending on the
particular application.
Absorbing Layers
Reference
1. J.S. Hesthaven and T. Warburton, Nodal Discontinuous Galerkin Methods —
Algorithms, Analysis, and Applications, Springer, 2008.
Initial Values
The Initial Values node adds the initial values for the Electric field and Magnetic field
variables that serve as an initial condition for the transient simulation.
INITIAL VALUES
Enter values or expressions for the initial values of the components of the Electric field
E (SI unit: V/m) and Magnetic field H (SI unit: A/m). The default values are 0 for all
vector components.
∂E
ε 0 ε r ------- – ∇ × H + σE = – J e
∂t
∂H
μ 0 μ r -------- + ∇ × E = – J m
∂t
Electric Field
The Electric Field boundary condition
n × E = n × E0
ELECTRIC FIELD
Enter values or expressions for the components of the Electric field E0 (SI unit: V/m).
n×E = 0
is a special case of the electric field boundary condition that sets the tangential
component of the electric field to zero. It is used for the modeling of a lossless metallic
surface, for example, a ground plane or as a symmetry type boundary condition.
It imposes symmetry for magnetic fields and antisymmetry for electric fields and
electric currents. It supports induced electric surface currents and thus any prescribed
or induced electric currents (volume, surface, or edge currents) flowing into a perfect
electric conductor boundary is automatically balanced by induced surface currents.
Magnetic Field
The Magnetic Field node adds a boundary condition for specifying the tangential
component of the magnetic field at the boundary:
n × H = n × H0
MAGNETIC FIELD
Enter values or expressions for the components of the Magnetic field H0
(SI unit: A/m).
n×H = 0
is a special case of the surface current density boundary condition that sets the
tangential component of the magnetic field and thus also the surface current density
–n × H = Js
n × ( H1 – H2 ) = Js
specifies a surface current density at both exterior and interior boundaries. The current
density is specified as a three-dimensional vector, but because it needs to flow along
the boundary surface, COMSOL Multiphysics projects it onto the boundary surface
and neglects its normal component. This makes it easier to specify the current density
and avoids unexpected results when a current density with a component normal to the
surface is given.
n × E = Z0 H
IMPEDANCE
Enter the value or expression for the medium Impedance Z0 (SI unit: Ω). By default,
the Z0 uses the value of the vacuum’s impedance. Then select Isotropic, Diagonal,
Symmetric, or Anisotropic based on the material characteristics and enter values or
expressions in the field or matrix.
n × E = E0
n × H = H0
specifies the tangential component of both electric and magnetic fields. This boundary
condition is available when Advanced Physics Options is checked from the Show menu
in the Model Builder toolbar.
BOUNDARY FLUX/SOURCE
Enter values or expressions for the components of the tangential Electric field E0
(SI unit: V/m) and the tangential Magnetic field H0 (SI unit: A/m).
Background Field
The Background Field feature triggers the scattered field formulation, where the
dependent variable is the relative field. The same wave equations are used as in the full
field formulation, but the total field that enters the equations are written as the sum of
the relative field and the background field, E = Erelative + Ebackground, and it is the
dependent variable Erelative that is solved for. When the background field is a solution
of the wave equation, the relative field is the scattered field.
SETTINGS
Select a Background wave type — User defined (the default), or Modulated Gaussian
pulse.
User Defined
Enter the component expressions for the Background electric field Eb (SI unit: V/m)
and Background magnetic field Hb (SI unit: A/m). The entered expressions must be
differentiable in time domain since the derivative of the background field is used in the
governing equations.
For a modulated Gaussian pulse propagating in the positive x-direction, the electric
field is expressed as
x + d offset 2
t – μ – ------------------------- -
vp
- sin 2πf 0 t – ------
1 - x
E ( x, t ) = ------------- exp – ---------------------------------------------------
τ 2π 2
v
2τ p
where τ is the pulse duration, defined as 1/2f0, μ is a time delay set to 2/f0, and vp is
the phase velocity. The time delay μ is used to excite a modulated Gaussian pulse whose
initial magnitude is very small when it is launched and gradually increases as it
propagates.
Far-Field Domain
To set up a far-field calculation, add a Far-Field Domain node and specify the far-field
domains in its Settings window. Use Far-Field Calculation subnodes (one is added by
default) to specify all other settings needed to define the far-field calculation. By
default, all of the domains are selected. The selection can be modified. In that case,
select only a homogeneous domain or domain group that is outside of all radiating and
scattering objects and which has the material settings of the far-field medium.
Far-Field Calculation
A Far-Field Calculation subnode is added by default to the Far-Field Domain node and
is used to select boundaries corresponding to a single closed surface surrounding all
radiating and scattering objects. By default, all exterior boundaries of the Far-Field
Domain are selected. Symmetry reduction of the geometry makes it relevant to select
boundaries defining a nonclosed surface. Also use this feature to indicate symmetry
planes and symmetry cuts applied to the geometry, and whether the selected
FAR-FIELD CALCULATION
Enter a Far-field variable name. The default is Efar.
Select as needed the Symmetry in the x=0 plane, Symmetry in the y=0 plane, or
Symmetry in the z=0 plane check boxes to use it your model when calculating the
far-field variable. The symmetry planes have to coincide with one of the Cartesian
coordinate planes.
When a check box is selected, also choose the type of symmetry to use from the
Symmetry type list that appears — Symmetry in E (PMC) or Symmetry in H (PEC). The
selection should match the boundary condition used for the symmetry boundary.
Using these settings, include the parts of the geometry that are not in the model for
symmetry reasons in the far-field analysis.
From the Boundary relative to domain list, select Inside or Outside (the default) to define
if the selected boundaries are defining the inside or outside of the far-field domain (that
is, whether facing away from infinity or toward infinity).
A Time to Frequency FFT study step must be added after the Time Dependent
study step to generate the necessary frequency-domain data, used in the
far-field analysis.
The physics interface can be used efficiently for unidirectional and bidirectional
propagation of electromagnetic beams. However, for optical scattering phenomena,
where the field is scattered into many different directions, the Electromagnetic Waves,
Frequency Domain interface is better suited.
With this physics interface the electric field is factored into a product of a slowly
varying envelope function (slowly on the scale of a wavelength) and a rapidly varying
phase function. The phase function is a priori prescribed, so the physics interface solves
the time-harmonic wave equation for the slowly varying envelope function.
The physics interface supports the study types Frequency domain, Wavelength
Domain, Eigenfrequency, and Boundary Mode Analysis. The frequency and
wavelength domain study types are used for source driven simulations for a single
frequency or wavelength or a sequence of frequencies or wavelengths. The
Eigenfrequency study type is used to find resonance frequencies and their associated
eigenmodes in cavity problems.
When this physics interface is added, these default nodes are also added to the Model
Builder—Wave Equation, Beam Envelopes, Perfect Electric Conductor, and Initial Values.
Then, from the Physics toolbar, add other nodes that implement, for example,
boundary conditions. You can also right-click Electromagnetic Waves, Beam Envelopes
to select physics features from the context menu.
SETTINGS
The Label is the default physics interface name.
The Name is used primarily as a scope prefix for variables defined by the physics
interface. Refer to such physics interface variables in expressions using the pattern
<name>.<variable_name>. In order to distinguish between variables belonging to
different physics interfaces, the name string must be unique. Only letters, numbers, and
underscores (_) are permitted in the Name field. The first character must be a letter.
COMPONENTS
This section is available for 2D and 2D axisymmetric models.
Select the Electric field components solved for—Three-component vector (the default),
Out-of-plane vector, or In-plane vector. Select:
• Three-component vector to solve using a full three-component vector for the electric
field envelope(s) E1 (and E2).
• Out-of-plane vector to solve for the electric field envelope vector component
perpendicular to the modeling plane, assuming that there is no electric field in the
plane.
• In-plane vector to solve for the electric field envelope vector components in the
modeling plane assuming that there is no electric field perpendicular to the plane.
WAVE VECTORS
Select the Number of directions—Bidirectional (the default) or Unidirectional.
Select the Type of phase specification—Wave vector (the default) or User defined.
In the tables, if Wave vector is selected for Type of phase specification, enter values or
expressions for the Wave vector, first wave k1 (SI unit: rad/m) and, if Bidirectional is
selected, for Wave vector, second wave k2 (SI unit: rad/m).
If User defined is selected for Type of phase specification, enter an expression for Phase,
first wave φ1 (SI unit: rad) and, if Bidirectional is selected, for Phase, second wave φ2
(SI unit: rad).
E ( r ) = E 1 ( r ) exp ( – jφ 1 ( r ) ) ,
where E1 is the electric field envelope that is solved for and exp(-jφ1(r)) is the
prescribed rapidly varying phase function. When Wave vector is selected for Type of
phase specification, the phase is defined as
φ1 ( r ) = k1 ⋅ r .
Notice that the wave vector k1 is assumed to be the same for all domains selected for
the physics interface. This also means that the phase will satisfy the condition of being
continuous everywhere. If the wave is assumed to bend or there are different materials
in the domains, the phase approximation above is not good and it is better to select a
The solution for the electric field envelope E1 is as exact as the solution for the total
electric field E, as is done for The Electromagnetic Waves, Frequency Domain
Interface. The advantage is that the mesh only need to resolve the spatial variation of
the field envelope E1 and not the rapid variation of the phase factor. On the other
hand, for problems involving reflections and scattering there is a rapid spatial variation
also for the field envelope. Then there is a no advantage of using the Unidirectional
formulation.
E ( r ) = E 1 ( r ) exp ( – jφ 1 ( r ) ) + E 2 ( r ) exp ( – jφ 2 ( r ) ) ,
where E2 and exp(-jφ2(r)) are the electric field envelope and the prescribed phase
function for the second wave. When specifying User defined phases, φ1 and φ2, each
phase should be continuous across the boundaries.
The Bidirectional formulation is good to use when there are boundaries reflecting the
wave in another direction than that of the incident wave. The direction for the reflected
beam is typically in the opposite direction to the incident beam. The boundary
conditions at these internal and/or external boundaries couple the electric field
envelopes E1 and E2.
Notice, however, that there is no coupling between E1 and E2 within domains, unless
weak expressions are explicitly added to the domains in the Model Builder. For more
information about how to add weak domain expressions, see Common Physics
Interface and Feature Settings and Nodes.
In the COMSOL Multiphysics Reference Manual see Table 2-3 for links
to common sections and Table 2-4 to common feature nodes. You can
also search for information: press F1 to open the Help window or Ctrl+F1
to open the Documentation window.
The default value for k2 represents the wave vector for a plane wave
reflected from a plane normal to the x-direction. Thus, the x-component
is negated, whereas the other components are the same as for wave vector
of the incident wave.
The default value for the User defined phase for the second wave, φ2,
represents a wave propagating in the opposite direction to the first wave.
For 2D axisymmetry, the default value for k1 (or φ1) represents a wave
vector pointing in the z-direction, whereas k2 (or φ2) represents a wave
propagating in the opposite direction to the first wave.
For an example using the User defined Type of phase specification, see
Gaussian Beam Incident at the Brewster Angle: Application Library
path Wave_Optics_Module/Optical_Scattering/brewster_interface.
The default values for the wave vectors are the gradients of the
corresponding phases defined in the Wave Vector settings. These values
will be correct for most cases. However, they will be wrong for Perfectly
Matched Layer domains. There, it is better to explicitly specify the wave
vector. For example, if the wave solution is expected to approximate a
plane wave in vacuum (or air), it would be better to enter the vacuum
wave number ewbe.k0 in the appropriate component field.
For this physics interface, the S-parameters are subject to Touchstone file export. Click
Browse to locate the file, or enter a filename and path. Select an Parameter format (value
pairs) — Magnitude angle, Magnitude (dB) angle, or Real imaginary.
Select an option from the If file exists list — Overwrite (the default) or Create new.
PHYSICS-CONTROLLED MESH
Select the Enable check box to use a physics-controlled mesh for the electromagnetic
problem. When Enable is selected, choose the Mesh type — Swept mesh (default for
3D), Mapped mesh (default for 2D), Tetrahedral mesh (3D), and Triangular mesh (2D).
When a structured Mesh type (either Swept mesh in 3D or Mapped mesh in 2D) is
selected, enter values for Number of transverse mesh elements (default is 10) and
Number of longitudinal mesh elements (default is 10). The entered Number of transverse
mesh elements will be distributed along the longest side of the input boundary. A
boundary is identified as an input boundary if there is an active feature, like a Port, a
Scattering Boundary Condition, etc, added to that boundary and the feature defines
an incident wave. The mesh will be denser in domains where the refractive index is
larger. Similarly, the entered Number of longitudinal mesh elements will be distributed
along propagation direction. Also here, the mesh will be denser in domains where the
refractive index is larger.
If no input features are defined, for instance for an eigenfrequency simulation, the
longitudinal direction is assumed to be the longest direction of the geometry and the
transverse plane is orthogonal to the longitudinal direction.
For an example using the Physics-controlled mesh with a Swept mesh, see
Directional Coupler: Application Library path
Wave_Optics_Module/Waveguides_and_Couplers/directional_coupler.
DEPENDENT VARIABLES
The dependent variables (field variables) are for the:
• Electric field envelope, first wave E1 and its components (in the Electric field envelope
components, first wave fields).
• Electric field envelope, second wave E2 and its components (in the Electric field
envelope components, second wave fields). The second wave is applicable if the Wave
Vectors are bidirectional.
The name can be changed but the names of fields and dependent variables must be
unique within a model.
DISCRETIZATION
DOMAIN
• Initial Values
• Polarization
• Wave Equation, Beam Envelopes
BOUNDARY CONDITIONS
With no surface currents present, the following boundary conditions for the electric
and magnetic fields need to be fulfilled
n II × ( E I – E II ) = 0
,
n II × ( H I – H II ) = 0
where the roman numerals denote the fields and normals on the two sides of the
boundary.
For the unidirectional formulation, the electric field is given by the product of the
electric field envelope E1 and the phase function (see Wave Vectors). Because E1 is
being solved for and the phase function is continuous across boundaries, the tangential
component of the electric field is always continuous, and thus the first condition is
automatically fulfilled. The second condition is equivalent to the natural boundary
condition for the unidirectional formulation
–1 –1
– n × [ ( μ r ∇ × E ) I – ( μ r ∇ × E ) II ] = n × jωμ 0 ( H I – H II ) = 0
For the bidirectional formulation the transverse electric and magnetic field envelopes
are not necessarily continuous across boundaries. Thus, the continuity of the transverse
electric and magnetic fields are enforced using weak expressions and constraints
applied on the boundary.
These features are also available and described for The Electromagnetic Waves,
Frequency Domain Interface:
E ( r ) = E 1 ( r ) exp ( – jϕ 1 ( r ) ) ,
for Wave Vectors set to unidirectional. Inserting this electric field formulation into the
Maxwell’s equations results in the following wave equation for the envelope function
2
( ∇ – jk 1 ) × ( ( ∇ – jk 1 ) × E 1 ) – k E 1 = 0 ,
where
k 1 = ∇ϕ 1 .
k = k0 n ,
where n is the refractive index and the wave number of free space k0 is defined as
ω
k 0 = ω ε 0 μ 0 = ----- .
c0
When Wave Vectors are set to bidirectional, the electric field is defined as the sum of
two fields
E ( r ) = E 1 ( r ) exp ( – jϕ 1 ( r ) ) + E 2 ( r ) exp ( – jϕ 2 ( r ) )
2
( ∇ – jk 2 ) × ( ( ∇ – jk 2 ) × E 2 ) – k E 2 = 0 ,
where
k 2 = ∇ϕ 2 .
ω
Q fact = --------- .
2δ
The settings for the Wave Equation, Beam Envelopes feature node are the
same as Wave Equation, Electric.
Polarization
The Polarization node adds an externally generated polarization Pi, which contributes
to the total polarization
P = ε 0 ( ε r – 1 )E + Pi .
i
P i = P 1i + P 2i .
Initial Values
The Initial Values node adds initial values for the electric field envelopes for the first and
second waves, which can serve as an initial guess for a nonlinear solver. Add more Initial
Values nodes from the Physics toolbar.
INITIAL VALUES
Enter values or expressions for the initial values of the components of the Electric field
envelope, first wave E1 and Electric field envelope, second wave E2 (SI unit: V/m). The
default values are 0 V/m. The second wave is applicable if the Wave Vectors are
bidirectional.
Electric Field
The Electric Field boundary condition
n × E = n × E0
specifies the tangential component of the electric field. It should in general not be used
to excite a model. Consider using the Port or Scattering Boundary Condition instead.
It is provided mainly for completeness and for advanced users who can recognize the
special modeling situations when it is appropriate to use. The commonly used special
case of zero tangential electric field is described in the Perfect Electric Conductor
section.
ELECTRIC FIELD
Enter the value or expression for the components of the Electric field, first wave E01
(SI unit: V/m). When Wave Vectors is set to bidirectional, also set the Electric field,
second wave E02.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
CONSTRAINT SETTINGS
To display this section, click the Show button ( ) and select Advanced Physics Options.
Magnetic Field
The Magnetic Field node adds a boundary condition for specifying the tangential
component of the magnetic field at the boundary:
n × H = n × H0
MAGNETIC FIELD
Enter the value or expression for the components of the Magnetic field, first wave H01
(SI unit: A/m). When Wave Vectors are set to bidirectional, also set the Magnetic field,
second wave H02.
Select an Incident field to specify whether the input wave is specified by the Electric field
or the Magnetic field.
Specify the Incident electric field envelope E0 (SI unit: V/m) or Incident magnetic field
envelope H0 (SI unit: A/m), depending on the Incident field selected. Notice that you
only specify the envelope factor of the incident electric or magnetic field. The envelope
function is internally multiplied by the phase function, as specified in the Wave Vectors
settings, to form the complete incident electric or magnetic fields.
When Wave Vectors is set to bidirectional, if no scattered field is expected, select the
No scattered field check box. This prevents COMSOL from returning spurious
solutions that otherwise could appear between boundaries with unconstrained
scattered fields.
The field E0 is the incident plane wave that travels in the direction k. The boundary
condition is transparent for incoming (but not outgoing) plane waves with any angle
of incidence. When Wave Vectors are set to unidirectional, the direction k is provided
automatically from the wave vector k1 specified for the physics interface. When Wave
Vectors are set to bidirectional, the user selects whether the direction k is provided
from the wave vector for the first wave k1 or the wave vector for the second wave k2.
• For cylindrical waves, specify around which cylinder axis the waves are cylindrical.
Do this by specifying one point at the cylinder axis and the axis direction.
• For spherical waves, specify the center of the sphere around which the wave is
spherical.
• When Wave Vectors are set to bidirectional, specify which wave the specified input
field is associated with.
– jk ( k ⋅ ( r – r ref ) )
E0 e ,
where rref is a reference point determined as the average point from the point selection
in the Reference Point subnode.
Select Incident field to specify whether the input wave is specified by the electric field
(Wave given by E field) or the magnetic field (Wave given by H field).
Specify the Incident electric field E0 (SI unit: V/m) or Incident magnetic field H0
(SI unit: A/m), depending on the setting of Incident field.
Select a Wave type for which the boundary is absorbing—Spherical wave, Cylindrical
wave, or Plane wave.
• For Cylindrical wave enter coordinates for the Source point ro (SI unit: m) and Source
axis direction raxis (dimensionless).
• For Spherical wave enter coordinates for the Source point ro (SI unit: m).
–n × H = Js
n × ( H1 – H2 ) = Js
specifies a surface current density at both exterior and interior boundaries. The current
density is specified as a three-dimensional vector, but because it needs to flow along
the boundary surface, COMSOL Multiphysics projects it onto the boundary surface
and neglects its normal component. This makes it easier to specify the current density
and avoids unexpected results when a current density with a component normal to the
surface is given.
∂D
∇ × H = J + -------
∂t
∂------
B-
∇×E = –
∂t
∂εE
∇ × H = σE + ----------
∂t
∂H
∇ × E = – μ --------
∂t
the two laws can be combined into a time harmonic equation for the electric field, or
a similar equation for the magnetic field
∇ × ( μ – 1 ∇ × E ) – ω 2 εE = 0
–1
∇ × ( ε ∇ × H ) – ω 2 μH = 0
The first of these, based on the electric field is used in The Electromagnetic Waves,
Frequency Domain Interface.
Using the relation εr = n2, where n is the refractive index, the equation can
alternatively be written
2 2
∇ × ( ∇ × E ) – k0 n E = 0 (3-1)
ω
k 0 = ω ε 0 μ 0 = -----
c0
When the equation is written using the refractive index, the assumption is that μr = 1
and σ = 0 and only the constitutive relations for linear materials are available. When
solving for the scattered field the same equations are used but E = Esc + Ei and Esc is
the dependent variable.
For The Electromagnetic Waves, Beam Envelopes Interface the electric field is written
as a product of an envelope function E1 and a rapidly varying phase factor with a
prescribed wave vector k1,
E ( r ) = E 1 ( r ) exp ( – jϕ 1 ( r ) ) .
When inserting this expression into Equation 3-1, the following wave equation for the
electric field envelope E1 is obtained
2 2
( ∇ – jk 1 ) × ( ( ∇ – jk 1 ) × E 1 ) – k 0 n E 1 = 0 , (3-2)
where
It is assumed that the envelope function E1 has a much slower spatial variation than
the exponential phase factor. Thus, the mesh can be much coarser when solving
Equation 3-2 than when solving Equation 3-1. Thereby it is possible do simulation on
domains that are much larger than the wavelength. Notice, however, that the
assumption of a slowly varying envelope function is never implemented in
Equation 3-2. Thus, the solution of Equation 3-2 is as exact as the solution of
Equation 3-1.
EIGENFREQUENCY ANALYSIS
When solving the frequency domain equation as an eigenfrequency problem the
eigenvalue is the complex eigenfrequency λ = jω + δ, where δ is the damping of the
solution. The Q-factor is given from the eigenvalue by the formula
ω
Q fact = ---------
2δ
The spatial parameter, α = δz + jβ = −λ, can have a real part and an imaginary part. The
propagation constant is equal to the imaginary part, and the real part, δz, represents
the damping along the propagation direction. When solving for all three electric field
components the allowed anisotropy of the optionally complex relative permittivity and
relative permeability is limited to:
PROPAGATING WAVES IN 2D
In 2D, different polarizations can be chosen by selecting to solve for a subset of the
3D vector components. When selecting all three components, the 3D equation applies
with the addition that out-of-plane spatial derivatives are evaluated for the prescribed
out-of-plane wave vector dependence of the electric field.
In 2D, the electric field varies with the out-of-plane wave number kz as (this
functionality is not available for The Electromagnetic Waves, Beam Envelopes
Interface)
˜
E ( x, y, z ) = E ( x, y ) exp ( – ik z z ) .
˜ 2 ˜
( ∇ – ik z z ) × [ μ r–1 ( ∇ – ik z z ) × E ] – k 0 ε rc E = 0 ,
Similarly, in 2D axisymmetry, the electric field varies with the azimuthal mode number
m as
˜
E ( r, ϕ, z ) = E ( r, z ) exp ( – imϕ )
∇ – im ˜ ˜
----- ϕ × μ r– 1 ∇ – i m
----- ϕ × E – k 20 ε rc E = 0 ,
r r
In-plane TM Waves
The TM waves polarization has only one magnetic field component in the z direction,
and the electric field lies in the modeling plane. Thus the time-harmonic fields can be
obtained by solving for the in-plane electric field components only. The equation is
formally the same as in 3D, the only difference being that the out-of-plane electric field
component is zero everywhere and that out-of-plane spatial derivatives are evaluated
for the prescribed out-of-plane wave vector dependence of the electric field.
In-plane TE Waves
As the field propagates in the modeling xy-plane a TE wave has only one nonzero
electric field component, namely in the z direction. The magnetic field lies in the
modeling plane. Thus the time-harmonic fields can be simplified to a scalar equation
for Ez,
2
– ∇ ⋅ ( μ̃ r ∇E z ) – ε rzz k 0 E z = 0
where
T
μr
μ̃ r = -------------------
det ( μ r )
Axisymmetric TM Waves
A TM wave has a magnetic field with only a ϕ component and thus an electric field
with components in the rz-plane only. The equation is formally the same as in 3D, the
only difference being that the ϕ component is zero everywhere and that spatial
derivatives with respect to ϕ are evaluated for the prescribed azimuthal mode number
dependence of the electric field.
Axisymmetric TE Waves
A TE wave has only an electric field component in the ϕ direction, and the magnetic
field lies in the modeling plane. Given these constraints, the 3D equation can be
simplified to a scalar equation for E ϕ . To write the fields in this form, it is also required
that εr and μr are nondiagonal only in the rz-plane. μr denotes a 2-by-2 tensor, and
ε rϕϕ and σ ϕϕ are the relative permittivity and conductivity in the ϕ direction.
Electric Losses
The frequency domain equations allow for several ways of introducing electric losses.
Finite conductivity results in a complex permittivity,
σ
ε c = ε – j ----
ω
ε c = ε 0 ( ε' – jε'' )
where ε' is the real part of εr, and all losses (dielectric and conduction losses) are given
by ε''. The dielectric loss model can also single out the losses from finite conductivity
(so that ε'' only represents dielectric losses) resulting in:
ε c = ε 0 ε' ( 1 – j tan δ )
For the physics interfaces in the Wave Optics Module, the refractive index is the default
electric displacement field model . In materials where μr is 1, the relation between the
complex refractive index
n = n – jκ
2
ε rc = n
that is
2 2
ε' r = n – κ
ε'' r = 2nκ
2 1 2 2
n = --- ( ε' r + ε' r + ε'' r )
2
2 1 2 2
κ = --- ( – ε' r + ε' r + ε'' r )
2
In the physics and optics literature, the time harmonic form is often written with a
minus sign (and “i” instead of “j”):
– i ωt
E ( x, y, z, t ) = E ( x, y, z )e
Magnetic Losses
The frequency domain equations allow for magnetic losses to be introduced as a
complex relative permeability.
μ r = ( μ' – jμ'' )
The complex relative permeability can be combined with any electric loss model except
refractive index.
∂A ∂ ∂A ∇ μ – 1
μ0 σ + μ0 ε + × ( r ∇ × A) = 0
∂t ∂t ∂t
Using the relation εr = n2, where n is the refractive index, the equations can
alternatively be written
2 ∂A
μ 0 ε 0 ∂ n + ∇ × (∇ × A) = 0
∂t ∂t
WAVES IN 2D
In 2D, different polarizations can be chosen by selecting to solve for a subset of the
3D vector components. When selecting all three components, the 3D equation applies
with the addition that out-of-plane spatial derivatives are set to zero.
In-plane TE Waves
As the field propagates in the modeling xy-plane a TE wave has only one nonzero
vector potential component, namely in the z direction. The magnetic field lies in the
modeling plane. Thus the equation in the time domain can be simplified to a scalar
equation for Az:
∂A z ∂ ∂A z
+ μ 0 ε 0 ε r
–1
μ0 σ + ∇ ⋅ ( μ r ( ∇A z ) ) = 0
∂t ∂t ∂t
Using the relation εr = n2, where n is the refractive index, the equation can
alternatively be written
∂ n 2 ∂A z
μ0 ε 0 + ∇ ⋅ ( ∇A z ) = 0
∂t ∂t
When using the refractive index, the assumption is that μr = 1 and σ = 0 and only the
constitutive relations for linear materials can be used.
Axisymmetric TM Waves
TM waves have a magnetic field with only a ϕ component and thus an electric field
and a magnetic vector potential with components in the rz-plane only. The equation
is formally the same as in 3D, the only difference being that the ϕ component is zero
everywhere and that spatial derivatives with respect to ϕ are set to zero.
Axisymmetric TE Waves
A TE wave has only a vector potential component in the ϕ direction, and the magnetic
field lies in the modeling plane. Given these constraints, the 3D equation can be
simplified to a scalar equation for A ϕ . To write the fields in this form, it is also required
Vector Elements
Whenever solving for more than a single vector component, it is not possible to use
Lagrange elements for electromagnetic wave modeling. The reason is that they force
the fields to be continuous everywhere. This implies that the physics interface
conditions, which specify that the normal components of the electric and magnetic
fields are discontinuous across interior boundaries between media with different
permittivity and permeability, cannot be fulfilled. To overcome this problem, the
Electromagnetic Waves, Frequency Domain interface uses vector elements, which do
not have this limitation.
The solution obtained when using vector elements also better fulfills the divergence
conditions ∇ · D = 0 and ∇ · B = 0 than when using Lagrange elements.
Eigenfrequency Calculations
When making eigenfrequency calculations, there are a few important things to note:
Using the default parameters for the eigenfrequency study, it might find a large
number of false eigenfrequencies, which are almost zero. This is a known consequence
of using vector elements. To avoid these eigenfrequencies, change the parameters for
the eigenvalue solver in the Study Settings. Adjust the settings so that the solver
searches for eigenfrequencies closer to the lowest eigenfrequency than to zero.
w0 2 2
ρ - ρ -
E b ( x, y, z ) = E bg0 ------------ exp – -------------- – jkz – jk --------------- + jη(z) ,
w( z) w (z)
2 2R ( z )
where w0 is the beam radius, p0 is the focal plane on the z-axis, Ebg0 is the background
electric field amplitude and the spot radius for different positions along the
propagation axis is given by
z – p0 2
w ( z ) = w 0 1 + --------------- .
z0
z0 2
R ( z ) = ( z – p 0 ) 1 + ---------------
z – p0
defines the radius of curvature for the phase of the field and the so called Gouy phase
shift is given by
z–p
η ( z ) = atan --------------0- .
z0
The equations above are expressed using the Rayleigh range z0 and the transverse
coordinate ρ, defined by
2
k0 w0 2 2 2
z 0 = --------------, ρ = x + y .
2
Note that the time-harmonic ansatz in COMSOL is ejωt and with this convention, the
beam above propagates in the +z-direction. The equations are modified accordingly
for beams propagating along the other coordinate axes.
The background field for a Gaussian beam is defined in a similar way for 2D
components. In the particular case where the beam propagates along the x-axis, the
background field is defined as
2 2
w0 y η(x)
y - -----------
E b ( x, y, z ) = E bg0 ------------ - – jkx – jk ---------------
exp – -------------- +j .
w(x ) w (x)
2 2R ( x ) 2
For a beam propagating along the y-axis, the coordinates x and y are interchanged.
To circumvent the problem that the paraxial approximation formula is not a solution
to the Helmholtz equation, a plane wave expansion can be used to approximate a
Gaussian beam. Since each plane wave is a solution to Helmholtz equation, also the
expansion is a solution to Helmholtz equation.
The plane wave expansion approximates the Gaussian distribution in the focal plane
x 2 + y 2
E b, Gauss ( r ) = E 0 exp – -----------------
2
- e =
w0
L M 1
where the beam is assumed to be propagating in the z-direction, the focal plane is
spanned by the x and y coordinates, e is the unit magnitude transverse polarization in
the focal plane, l and m denote the indices for the wave vectors, the index n accounts
for the two polarizations per wave vector klm, almn is the amplitude, un(klm) is the
unit magnitude polarization, and r is the position vector.
Multiplying with the conjugate of the exponential factor above and the polarization
factor un(klm) and applying a surface integral over the entire focal plane allows us to
extract the amplitudes as
2 2 2
E 0 w 0 ( e ⋅ u n ( k lm ) ) k t, lm w 0
a lmn = -------------------------------------------------- exp – --------------------- ,
4π 4
• The Equations
• In-plane E Field or In-plane H Field
• Fluxes as Dirichlet Boundary Conditions
• Absorbing Layers
The Equations
Maxwell’s equations are a set of equations, written in differential or integral form,
stating the relationships between the fundamental electromagnetic quantities. These
quantities are the:
For general time-varying fields, the differential form of Maxwell’s equations can be
written as
∂D
∇ × H = J + -------
∂t
∂B
∇ × E = – ------- (3-3)
∂t
∇⋅D = ρ
∇⋅B = 0
The first two equations are also called Maxwell-Ampere’s law and Faraday’s law,
respectively. Equation three and four are two forms of Gauss’ law, the electric and
magnetic form, respectively.
D = ε0 E + P
B = μ0 ( H + M ) (3-4)
J = σE
Here ε0 is the permittivity of a vacuum, μ0 is the permeability of a vacuum, and σ the
electric conductivity of the medium. In the SI system, the permeability of a vacuum is
chosen to be 4π·10−7 H/m. The velocity of an electromagnetic wave in a vacuum is
given as c0 and the permittivity of a vacuum is derived from the relation
1 – 12 1 –9
ε 0 = ----------
2
= 8.854 ⋅ 10 F/m ≈ --------- ⋅ 10 F/m
c0 μ0 36π
The electric polarization vector P describes how the material is polarized when an
electric field E is present. It can be interpreted as the volume density of electric dipole
moments. P is generally a function of E. Some materials might have a nonzero P also
when there is no electric field present.
The magnetization vector M similarly describes how the material is magnetized when
a magnetic field H is present. It can be interpreted as the volume density of magnetic
dipole moments. M is generally a function of H. Permanent magnets, for example,
have a nonzero M also when there is no magnetic field present.
To get a wave equation for the E field, for example, take the curl of the second
equation in Equation 3-3 (previously divided by μ0), and insert it into the time
derivative of the first row in Equation 3-3
2 2
∂M ∂E ∂ E- ∂ P
– ∇ × ----- ∇ × E + -------- = σ ------- + ε 0 ---------
1
+ ---------
-
μ0 ∂t ∂t ∂t 2 ∂t 2
this is referred as curl-curl formulation in the literature (second order time derivatives
and second order space derivatives).
LINEAR MATERIALS
In the simplest case linear materials, the polarization is directly proportional to the
electric field, that is
∂P ⁄ ∂E = ε 0 χ e and P = ε 0 χ e E
∂M ⁄ ∂H = χ m and M = χ m H
As a consequence, for linear materials, the constitutive relations in Equation 3-4 can
be written as
D = ε 0 E + P = ε 0 ( 1 + χ e )E = ε 0 ε r E
B = μ 0 ( H + M ) = μ 0 ( 1 + χ m )H = μ 0 μ r H
Here, ε = ε0εr and μ = μ0μr are the permittivity and permeability of the material. The
relative permittivity εr and the relative permeability μr are usually scalar properties but
these can be second-rank symmetric (Hermitian) tensors for a general anisotropic
material.
∂E
∇ × H = σE + ε 0 ε r -------
∂t
(3-5)
∂-------
H-
∇ × E = –μ0 μr
∂t
∂u
da + ∇ ⋅ Γ(u) = f
∂t
Maxwell’s equations in 3D
∂E
ε 0 ε r ------- – ∇ × H = – σE
∂t
∂H
μ 0 μ r -------- + ∇ × E = 0
∂t
∂E
d E ------- + ∇ ⋅ Γ E ( H ) = f
∂t
∂-------
H
dH - + ∇ ⋅ ΓH ( E ) = 0
∂t
d E = ε 0 ε r and d H = μ 0 μ r
0 h3 –h2 0 e3 –e2
ΓE ( H ) = – –h3 0 h1 and Γ H ( E ) = – e 3 0 e1
h2 –h1 0 e2 –e1 0
When using SI units (or other) for the electromagnetic fields and material
properties, the Lax-Friedrichs flux parameter is not dimensionless and
must have units of τE = 1/(2Z) for Ampere’s law and τH = Z/2 for
Faraday’s law, where Z is the impedance of the medium.
TM WAVES IN 2D
For TM waves in 2D, solve for an in-plane electric field vector and one out-of-plane
variable for the magnetic field. Maxwell’s equations then read
∂E
ε 0 ε r ------- + ∇ ⋅ Γ E ( H ) = – σ ⋅ E
∂t
(3-7)
∂H
μ 0 μ r -------- + ∇ ⋅ Γ H ( E ) = 0
∂t
0 –h3
ΓE ( H ) = and Γ H ( E ) = e 2 – e 1 (3-8)
h3 0
The default Lax-Friedrichs flux parameters are τE = 1/(2Z) for Ampere law, and the
scalar τH = Z/2 for Faraday’s law, where Z is the impedance of a vacuum.
TE WAVES IN 2D
For TE waves in 2D, solve for an in-plane magnetic field vector and one out-of-plane
variable for the electric field. Maxwell’s equations then read
∂E
ε 0 ε r ------- + ∇ ⋅ Γ E ( H ) = – σE
∂t
(3-9)
∂H
μ 0 μ r -------- + ∇ ⋅ Γ H ( E ) = 0
∂t
0 e3
Γ E ( H ) = – h 2 h 1 and Γ H ( E ) = (3-10)
–e3 0
The default Lax-Friedrichs flux parameters are τE = 1/(2Z) for Ampere law, and two
scalar τH = Z/2 for Faraday’s law, where Z is the impedance of a vacuum.
∂E
ε 0 ε r ------- + ∇ ⋅ Γ E ( H ) = – σE
∂t
∂H
μ 0 μ r -------- + ∇ ⋅ Γ H ( E ) = 0
∂t
0 –h3 h2 0 e3 –e2
ΓE ( H ) = h3 0 – h 1 and Γ H ( E ) = – e 3 0 e1
–h2 h1 0 e2 –e1 0
For Ampere’s law, the normal to the flux term on exterior boundaries reads
n ⋅ ΓE ( H ) = –n × H
n ⋅ ΓH ( E ) = n × E
which means that normal fluxes on external boundaries can only prescribe tangential
components for the fields.
BOUNDARY CONDITIONS
The boundary conditions for outer boundaries are computed from the normal fluxes
n · ΓH(E) and n · ΓE(H).
Absorbing Layers
The Electromagnetic Waves, Time Explicit Interface includes so-called absorbing
layers, also often referred to as sponge layers. The layers work by combining three
techniques: a scaling system, filtering, and simple nonreflecting conditions. For a
review of the method see, for example, Ref. 1.
The layers are set up by adding the Absorbing Layer under the Definitions node. This
adds a special scaled system. The scaling effectively slows down the propagating waves
and ensures that they hit the outer boundary in the normal direction. For the
Absorbing Layer domain selection, add an additional Wave Equations feature, mark
the Activate check box under the Filter Parameters section, and enter filter parameters.
Filtering attenuates and filters out high-frequency components of the wave. Finally, at
the outer boundary of the layer add a simple Scattering Boundary Condition
condition, which will work well to remove all remaining waves as normal incidence has
been ensured.
For more detailed information about the filter see the Filter Parameters
section under Wave Form PDE in the COMSOL Multiphysics Reference
Manual.
For more detailed on the Geometry and Scaling see the Infinite Elements,
Perfectly Matched Layers, and Absorbing Layers in the COMSOL
Multiphysics Reference Manual.
For the layers to work optimally the filter should not be too aggressive. Moreover, the
scaled coordinates in the layer domain should also vary smoothly. To inspect the scaled
system you can plot the coordinate variables x_absorb_ab1, y_absorb_ab1, and
z_absorb_ab1. Using the absorbing layers with the three combined techniques will
enable the reduction of spurious reflections by a factor between 100 and 1000
compared to the incident amplitude.
Reference
1. P.G. Petropoulos, L. Zhao, and A.C. Cangellaris, “A Reflectionless Sponge Layer
Absorbing Boundary Condition for the Solution of Maxwell’s Equations with
High-Order Staggered Finite Difference Schemes,” J. Comp. Phys., vol. 139, pp. 184–
208, 1998.
This chapter describes The Laser Heating Interface found under the Heat
Transfer>Electromagnetic Heating branch ( ) when adding a physics interface.
See The Heat Transfer Interfaces and The Joule Heating Interface in the COMSOL
Multiphysics Reference Manual for other Heat Transfer interface and feature node
settings.
165
The Laser Heating Interface
The Laser Heating interface ( ) is used to model electromagnetic heating for systems
and devices where the electric field amplitude varies slowly on a wavelength scale. This
multiphysics interface adds an Electromagnetic Waves, Beam Envelopes interface and
a Heat Transfer in Solids interface. The multiphysics couplings add the
electromagnetic losses from the electromagnetic waves as a heat source, and the
electromagnetic material properties can depend on the temperature. The modeling
approach is based on the assumption that the electromagnetic cycle time is short
compared to the thermal time scale.
The Heat Transfer in Solids interface provides features for modeling heat transfer by
conduction, convection, and radiation. A Heat Transfer in Solids model is active by
default on all domains. All functionality for including other domain types, such as a
fluid domain, is also available. The temperature equation defined in solid domains
corresponds to the differential form of Fourier's law that may contain additional
contributions like heat sources.
However, if physics interfaces are added one at a time, followed by the coupling
features, these modified settings are not automatically included.
For example, if single Electromagnetic Waves, Beam Envelopes and Heat Transfer in Solids
interfaces are added, COMSOL Multiphysics adds an empty Multiphysics node. You
can choose from the available coupling features, Electromagnetic Heat Source and
Temperature Coupling, but any modified settings are not included.
Coupling features are available from the context menu (right-click the
Multiphysics node) or from the Physics toolbar, Multiphysics menu.
Electromagnetic Heat The Domain Selection is the same as that of the participating
Source physics interfaces.
The Boundary Selection is the same as the exterior and interior
boundaries of the Domain Selection of the participating physics
interfaces.
The corresponding Electromagnetic Waves, Beam Envelopes and
Heat Transfer in Solids interfaces are preselected in the Coupled
Interfaces section (described in the COMSOL Multiphysics
Reference Manual).
Temperature Coupling The corresponding Electromagnetic Waves, Beam Envelopes and
Heat Transfer in Solids interfaces are preselected in the Coupled
Interfaces section (described in the COMSOL Multiphysics
Reference Manual).
A side effect of adding physics interfaces one at a time is that four study
types — Frequency-Stationary, Frequency-Transient, Sequential
Frequency-Stationary, and Sequential Frequency-Transient— are not
available for selection until after at least one coupling feature is added. In
this case, it is better to initially not add any study at all, then add the
coupling features to the Multiphysics node, and lastly, open the Add Study
window and add a study sequence below the Preset Studies for Selected
Multiphysics heading.
Use the online help in COMSOL Multiphysics to locate and search all the
documentation. All these links also work directly in COMSOL
Multiphysics when using the Help system.
Coupling Features
The Electromagnetic Heating and Temperature Coupling coupling feature nodes are
described for The Joule Heating Interface in the COMSOL Multiphysics Reference
Manual.
• The available physics features for The Electromagnetic Waves, Beam Envelopes
Interface are listed in the section Domain, Boundary, Edge, and Point Nodes for
the Electromagnetic Waves, Beam Envelopes Interface.
• The available physics features for The Heat Transfer Interfaces are listed in the
section Feature Nodes for the Heat Transfer in Solids Interface in the COMSOL
Multiphysics Reference Manual.
Glossary
171
Glossary of Terms
absorbing boundary A boundary that lets an electromagnetic wave propagate through
the boundary without reflections.
constitutive relation The relation between the D and E fields and between the B and
H fields. These relations depend on the material properties.
cutoff frequency The lowest frequency for which a given mode can propagate
through, for example, a waveguide or optical fiber.
eigenmode A possible propagating mode of, for example, a waveguide or optical fiber.
electric dipole Two equal and opposite charges +q and −q separated a short distance
d. The electric dipole moment is given by p = qd, where d is a vector going from −q
to +q.
magnetic dipole A small circular loop carrying a current. The magnetic dipole
moment is m = IAe, where I is the current carried by the loop, A its area, and e a unit
vector along the central axis of the loop.
perfect electric conductor (PEC) A material with high electrical conductivity, modeled
as a boundary where the electric field is zero.
surface current density Current density defined on the surface. The component
normal to the surface is zero. The unit is A/m.
vector element A finite element often used for electromagnetic vector fields. The
tangential component of the vector field at the mesh edges is used as a degree of
freedom. Also called Nedelec’s edge element or just edge element.
INDEX| 175
constitutive relations 157 electromagnetic energy theory 44
constitutive relations, theory 42 electromagnetic quantities 64
continuity, periodic boundaries and 30 electromagnetic sources, applying 28
curl-curl formulation 157 electromagnetic waves, beam envelopes
cutoff frequency 51 interface 129
cylindrical coordinates 26 electromagnetic waves, frequency do-
cylindrical waves 98, 141 main interface 68–69
theory 144
D Debye dispersion model 80
electromagnetic waves, time explicit in-
dielectric medium theory 48
terface 118
diffraction order (node) 93
theory 156
dispersive materials 45
electromagnetic waves, transient inter-
documentation 17
face 110
domain nodes
theory 144
electromagnetic waves, beam enve-
emailing COMSOL 19
lopes 134
exponential filter, for wave problems 121
electromagnetic waves, frequency do-
external current density (node) 82
main interface 75
electromagnetic waves, time explicit F far field variables 36
119 Faraday’s law 156
Drude-Lorentz dispersion model 80 far-field calculation (node) 83, 127
Drude-Lorentz polarization (node) 116 far-field calculations 53
far-field domain (node) 82, 127
E E (PMC) symmetry 34
far-field variables 34
edge current (node) 108
field continuity (node) 139
edge nodes
file, Touchstone 74, 133
electromagnetic waves, beam enve-
Floquet periodicity 30, 107
lopes 134
flux/source (node) 126
eigenfrequency analysis 57
free-space variables 77, 136
eigenfrequency calculations theory 153
frequency domain equation 144
eigenfrequency study 146
Frequency-Domain Modal Method 63
eigenmode analysis 50
eigenvalue (node) 58 G Gauss’ law 156
electric current density (node) 123
H H (PEC) symmetry 34
electric field (node) 96, 123, 138
hybrid-mode waves
electric losses theory 149
axisymmetric, frequency domain 149
electric point dipole (node) 108
axisymmetric, time domain 152
electric scalar potential 44
in-plane, frequency domain 148
electric susceptibility 158
in-plane, time domain 151
electrical conductivity 42
176 | I N D E X
perpendicular 147 mesh resolution 28
mode analysis 59, 146
I impedance boundary condition (node)
mode phase
101
for Port and Diffraction Order 89
inhomogeneous materials 45
modeling tips 24
initial values (node)
MPH-files 18
electromagnetic waves, beam enve-
lopes 138 N nonlinear materials 45
electromagnetic waves, frequency do- numeric modes 87
main interface 82
P pair nodes
electromagnetic waves, time explicit
electromagnetic waves, beam enve-
interface 122
lopes 134
electromagnetic waves, transient 116
PEC. see perfect electric conductor
in-plane TE waves theory
perfect conductors theory 48
frequency domain 148
perfect electric conductor (node) 124
time domain 152
boundaries 84
in-plane TM waves theory
perfect magnetic conductor (node) 85,
frequency domain 148
124
time domain 152
periodic boundary conditions 30
inports 88
periodic condition (node) 106
internet resources 16
periodic port reference point (node) 95
K knowledge base, COMSOL 19 permeability
anisotropic 146
L line current (out-of-plane) (node) 109
permeability of vacuum 42
linearization point 58
permittivity
listener ports 88
anisotropic 146
losses, electric 149
permittivity of vacuum 42
losses, magnetic 151
phasors theory 48
M magnetic current (node) 107
PMC. see perfect magnetic conductor
magnetic current density (node) 123
point nodes
magnetic field (node) 97, 124, 139
electromagnetic waves, beam enve-
magnetic losses theory 151
lopes 134
magnetic point dipole (node) 108
polarization (node) 84, 137
magnetic susceptibility 43, 158
polarization, 2D and 2D axisymmetry 26
matched boundary condition (node) 139
port (node) 86
material properties 45
port boundary conditions 56
materials 46
potentials theory 44
Maxwell’s equations 41
Poynting’s theorem 44
Maxwell-Ampere’s law 156
INDEX| 177
Q quality factor (Q-factor) 57, 146 symmetry planes, far-field calculations 34
178 | I N D E X
S-parameters 56
vector elements theory 153
INDEX| 179
180 | I N D E X