# Assessment of a PML Boundary Condition for Simulating an MRI Radio Frequency Coil

Authors: Yunsuo Duan; Tamer S. Ibrahim; Bradley S. Peterson; Feng Liu; Alayar Kangarlu
Journal: International Journal of Antennas and Propagation
 (2008)

---

## Abstract

Computational methods such as the finite difference time domain (FDTD)  play an important role in simulating radiofrequency (RF) coils used in magnetic resonance imaging (MRI). The choice of absorbing boundary conditions affects the final outcome of such studies. We have used FDTD to assess the Berenger's perfectly matched layer (PML) as an absorbing boundary condition for computation of the resonance patterns and electromagnetic fields of RF coils. We first experimentally constructed a high-pass birdcage head coil, measured its resonance pattern, and used it to acquire proton(1H) phantom MRI images. We then computed the resonance pattern and B1 field of the coil using FDTD with a PML as an absorbing boundary condition. We assessed the accuracy and efficiency of PML by adjusting the parameters of the PML and comparing the calculated results with measured ones. The optimal PML parameters that produce accurate (comparable to the experimental findings) FDTD calculations are then provided for the birdcage head coil operating at 127.72 MHz, the Larmor frequency of 1H at 3 Tesla (T).


---

## Body

## 1. Introduction

The dimensions and resolutions
of a discretized domain for calculating electromagnetic fields can be
restrained by the memory and computational capacity of the computer. As a
result, the computational domain must utilize efficient and accurate boundary
conditions in order to truncate outward waves and therefore simulate infinite
radiation boundary conditions of the computational domain. One of the most useable
and efficient absorbing boundary conditions that have been reported for the finite
difference time domain (FDTD) [1] is the perfect matched layer (PML) [2].Since its introduction to FDTD, PML has been extended by a number of
researchers. Among these PMLs, Berenger’s PML [2, 3], anisotropic PML (APML)
[4, 5], complex frequency shifted PML (CFS PML) [6], and higher-order
PML [7, 8] represent the most usual types of PML. The significant distinction
between Berenger’s PML and APML, which together are designated regular PML, is
that the former uses Maxwellian isotropic absorbing materials and splits each
component of electromagnetic fields into two components, whereas the latter uses
anisotropic absorbing materials and keeps the components of electromagnetic
fields unsplit. Thus, the implementation of Berenger’s PML requires more computer
RAM than does APML. CFS PML is more efficient in annihilating the evanescent
waves than are the regular PMLs [9], although it sacrifices good absorption of the
propagation waves at low frequencies [7]. Higher-order PML is a new attempt to
combine the advantages of both regular PML and CFS PML [7, 8].PMLs have been widely used
with FDTD to compute electromagnetic fields, including those of magnetic
resonance imaging (MRI) [10–13]. Thus far, however, the assessment
of the accuracy and efficiency of PML when applied to MRI radiofrequency (RF)
simulations has not been reported. In fact, the required grid sizes conflict with
the required resolutions of the calculated frequency when using FDTD to compute
the resonance pattern and electromagnetic fields of RF coils. On the one hand, the
small dimensions of copper strips (~30 μm thick) used to construct RF coils
require small grids matched to the thickness of the copper strips so as to
improve the accuracy of the simulations. On the other hand, the frequency
resolution of the calculated resonance patterns and the requisite computing
time are inversely proportional to the size of the grids. Thus, the grid sizes
should be sufficiently large to improve the frequency resolution of the
calculated resonance pattern of the coil and to reduce computing time. Moreover,
because the dimensions of PML together with the coil geometry are enormous relative
to the size of the copper strips, a vast number of Yee cells are needed to
simulate the coil fully. Satisfying these conflicting requirements requires both
a vast computer RAM and an enormous number of iterated time steps, which
produces an unacceptable accumulation of numerical phase errors as well as an unrealistic
computational time. A PC with 3.39 GHz dual processors and 2 GB RAM, for
example, requires about 116 hours to compute the resonance pattern of an adult
RF head coil using FDTD with a grid size
of 2 × 2 × 2 mm3 and frequency
resolution of 0.5 MHz. If the calculated resonance pattern does not match the
desired one (the coil is off tuned), this computation must be repeated using
modified coil parameters until the desired resonance pattern is achieved.Furthermore, the increasing
strength of the static magnetic fields (B0) in MRI applications has
reduced the RF wavelength inside human tissue close to or below the dimensions
of the subject to be scanned. For example, at 3 Tesla (T), the wavelength of an
RF field inside the human head is about 280 mm, which is close to the average
diameter of an adult RF head coil. Therefore, the geometric structure of RF coils
and the sizes of the biological load (head or body) are no longer negligible in
relation to the wavelength of the operational frequency when B0 is 3 T
or higher. Thus, conventional lumped-element methods such as transmission line and
circuit theories become inadequate for analyzing and designing RF coils. As an
alternative, the FDTD method, a full-wave numerical technique that is capable
of accounting for the geometry of a coil as well as the intricacies of the sample,
has emerged as the tool of choice for analyzing RF coils at high fields [10–13]. Consequently,
the need for a good computationally efficient boundary condition has become
more pressing.We present a systematic
evaluation of the effect of the parameters of Berenger’s PML on the accuracy
and efficiency of the calculated resonance pattern and the electromagnetic
fields of RF coils. Because the Larmor
frequency of1H at 3 T (127.72 MHz) is much lower than the
frequencies at which the CFS PML works well, we have omitted the use of CFS PML
for simulation of RF coils. Furthermore, because the additional RAM required for
Berenger’s PML compared to APML is not an excessive burden for modern computers,
and because implementation of Berenger’s PML is easier than that of APML, we selected
Berenger’s PML as the absorbing boundary condition for assessment of simulated
RF coils. This choice does not limit the
generality of our assessments for other PMLs, because the influence trends of the
parameters of other PMLs are similar to those of Bérenger’s [9].We first constructed an actual high-pass birdcage head coil [14], measured
its resonance pattern, and acquired proton images from a phantom. We then computationally
modeled the resonance pattern and B1 fields of the coil (i.e., the circularly
polarized component of the transverse magnetic field that is responsible for
excitation) using an FDTD code developed in-house that uses Berenger’s PML as
an absorbing boundary condition. To assess the accuracy and efficiency of PML
for this purpose, we investigated the criteria for optimizing PML parameters by
comparing the numerically calculated results with the experimentally measured ones
until the best agreement between them was achieved. Thus, we identified the optimal
PML parameters that achieved both high accuracy and at least moderate computational
efficiency for analyzing the essential performance characteristics of RF coils.


## 2. Material and Methods

Our method for
assessing PML differed from the approaches used in the design of RF coils. Typically,
those approaches first evaluate theB1 field of a modeled coil and then
optimize it by changing the geometry and lumped components of the coil’s design
before constructing it. Here, our goal was not to design a coil, but to assess the
performance of the PML in accurately and efficiently simulating an RF coil.
Therefore, we first designed, constructed, and tuned a birdcage head coil to
the proper operational mode. We measured the coil’s resonance pattern and used
it to acquire 1H images of a phantom in a 3 T MRI scanner. Next, this
actual coil was modeled using an FDTD scheme that used PML as an absorbing
boundary condition. We calculated the resonance pattern and B1 field
of the coil. We compared the calculated and measured results to adjust the
parameters of the PML and to assess the accuracy of the PML according to the criteria
defined in the following sections. We further assessed computational efficiency
of the PML by comparing the computational demands for each set of PML parameters,
while assuming that greater efficiency was associated with less computational time.


### 2.1. Experimental Measurements

Our experimental head coil was constructed based on the high-pass
birdcage design [14] for a typical adult head. The coil had 16 rungs connected
by 2 end-rings, 32 capacitors of 13.6 pF, a diameter of 292 mm, and a length of
356 mm. The rungs and rings were made of copper strips having a width of 8 mm and
a thickness of 30 um (see Figure 1).(a) Structure of a high-pass birdcage head coil
showing a schematic view of the coil and (b) coronal and sagittal planes of the
coil, where “l”: length of the coil; “d”:
diameter of the coil; “C”: capacitor; “A”
and “B”: RF feeding points; “O”:
center of the coil.
(a)(b)We obtained the resonance pattern of the coil by measuring
the S-parameters using an Agilent 4395A network/spectrum/impedance analyzer. The
coil worked in quadrature mode [14]
with a 90° hybrid coupler
connected between the analyzer and the coil’s excitation ports. We acquired 1H
images from a phantom loaded in the coil using a GE 3 T Signa scanner and a
gradient echo pulse sequence having a 90° flip angle, 500 milliseconds
repetition time (RT]), and 10 milliseconds echo time (ET).The phantom consisted of a spherical plastic container filled
with purified water and NaCl in order to approximately simulate the average
dielectric constant (70ε0) and conductivity (0.57 S/m) of the human
brain at 128 MHz [9, 15]. The corresponding wavelength inside the phantom was
280 mm, which was 95.6% of the diameter and 78.6% of the length of the actual coil.


### 2.2. Computational Model

In simulating
the RF coil, we first used FDTD with Berenger’s PML as the boundary condition
to model the actual coil loaded with the phantom. We then computed the
resonance pattern andB1 field of the actual coil using the modeled
coil. We finally evaluated performance of the PML by comparing the calculated
results with the experimental ones. Our FDTD model and computing procedures are
described below.


#### 2.2.1. FDTD Model

Our computer code for FDTD was developed in-house. The FDTD model
had 3D grids of 2 mm definition in each of thex-, y-, and z-directions. The
z-axis of the simulation was parallel to the longitudinal axis of the coil. Although
the actual coil and the phantom were modeled according to their specific
physical sizes, some of the modeled dimensions deviated from those of the real
coil because the dimensions of some of the real components were smaller than the
definition of our grid. For example, the copper strips were represented by one
voxel (2 mm) in thickness rather than by their actual thickness of 30 microns. Similarly,
the lumped capacitors were represented by one voxel and modeled by the
following equation [16]:
(1)Ixn+1/2=CΔxExi,j,kn+1−Exi,j,kn−12Δt,
where Ix was the current flowing
though the capacitor, C was the capacitance of the capacitor, and x was the
direction along which the capacitor aligned. Note that “x” was substituted with
“y” or “z” when the capacitor aligned along the respective directions.For stable time marching, the time step length,Δt must satisfy the Courant-Friedrichs-Lewy (CFL)
stability criterion [1]:
(2)Δt<cmax⁡1/Δx2+1/Δy2+1/Δz2,
where cmax⁡ is the maximum propagation velocity of the electromagnetic waves within the
simulated region. In this case, Δx=Δy=Δz=2⁢mm,
and cmax⁡, was the velocity
of light. The corresponding maximum  Δt was 3.789 picoseconds.
According to the Nyquist criterion for sampling, the associated frequency resolution
(Δf) was given as
(3)Δf=1(NΔt),
where N was the number of time steps (iterations).
For a frequency resolution of 1 MHz or less, N should be greater than 263922.
We selected N as 262144 (= 218) for convenience of calculating the fast
Fourier transformation (FFT).


#### 2.2.2. Computation of Resonance Pattern

We
computed the resonance pattern of the coil by applying the FFT to its calculated
magnetic field. In so doing, we first excited the coil with an excitation signal
that contained a sufficient number of frequencies such that the calculated
electromagnetic fields of the coil contained all frequencies within its bandwidth.
Typically, the resonance frequency of the highest mode of the birdcage
coil for1H at 3 T is <500 MHz. We, therefore, chose a Gaussian
function as the excitation signal having a −3 dB frequency of 500 MHz:
(4)s(t)=e−(t/2T0−3)2,
where T0=1/ω−3⁢dB=1/(109π).We then applied this excitation signal to a capacitor of the coil and
iteratively computed the electromagnetic field of the coil. We saved the
calculatedB1 field for every time step at sampling points on the
surface and at the center of the coil (see Figure 1(b)). We then applied FFT to
the saved B1 field to obtain the frequency spectra (i.e., the resonance
pattern) of the coil at the related sampling points.A correctly calculated resonance pattern should agree with the real coil’s
actual resonance pattern provided that the frequency deviations between the
calculated and measured patterns are below the resolution of the calculated
resonance pattern and therefore are considered negligible. Otherwise, we must adjust the parameters
of the PML or the parameters of the modeled coil and recompute the resonance
pattern.


#### 2.2.3. Computation of theB1 Field

After calculating the coil’s resonance
pattern, we then computed theB1 field of the coil at the Larmor
frequency of 1H at 3 T. To achieve quadrature excitation, we,
respectively, applied two sine waves centered at the Larmor frequency to feed two
capacitors 90° azimuthally apart from each other (points A and B in Figure 1(b)). The two sine wave excitations having a 90° difference in their
phases were presented as
(5a)s1(t)=(1−e−(ω0t/5)2)sin⁡(2πf0t),(5b)s2(t)=(1−e−(ω0t/5)2)sin⁡(2πf0t−π/2),
where f0 is the Larmor frequency (127.72 MHz for 1H
at 3 T). The introduction of the coefficient (1−e−(ω0t/5)2) in the excitation signal accelerated
stabilization of the calculated B1 field [10].


### 2.3. PML

#### 2.3.1. PML Parameters

When implementing Berenger’s PML [2, 3], layers
of artificial perfect electric conductor (PEC) were introduced at the boundaries
of the computational domain to absorb the waves traveling outwards (i.e., to satisfy
the radiation condition). The conductivity of each layer increased gradually
from inside (near the coil’s geometry) to outside (away from the coil’s
geometry). Thus, the outgoing waves attenuated to negligible levels as they
traveled through the PML, then reflected at the PEC walls located at the outer
boundaries of the computational domain, and then traveled through the PML back
to the geometry of the coil.To compute the fields in the presence of Berenger’s
PML, each component of the electromagnetic fields must be divided into two
subcomponents. For example, the electric field componentEx was
divided into Exy and Exz:
(6)Ex=Exy+Exz,
where
(7a)Exyn+1(i+12,j,k)=e−(σy/ε0)ΔtExyn(i+12,j,k)+1σyΔy(1−e−(σy/ε0)Δt)×{[Hzxn+1/2(i+12,j+12,k)−Hzxn+1/2(i+12,j−12,k)]⁢−[Hzyn+1/2(i+12,j,k+12)−Hzyn+1/2(i+12,j,k−12)]},(7b)Exzn+1(i+12,j,k)=e−(σz/ε0)ΔtExzn(i+12,j,k)+1σzΔz(1−e−(σz/ε0)Δt)×{[Hyxn+1/2(i+12,j+12,k)−Hyxn+1/2(i+12,j−12,k)]⁢−[Hyzn+1/2(i+12,j,k+12)−Hyzn+1/2(i+12,j,k−12)]}.
Similarly, the magnetic field component Hx was divided
into Hxy and Hxz:
(8)Hx=Hxy+Hxz,
where
(9a)Hxyn+1/2(i,j+12,k+12)=e−(σy*/μ0)ΔtHxyn−1/2(i,j+12,k+12)+1σy*Δy(1−e−(σy*/μ0)Δt)×{[Ezxn(i,j+1,k+12)−Ezxn(i,j,k+12)]⁢−[Ezyn(i,j+1,k+12)−Ezyn(i,j,k+12)]},(9b)Hxzn+1/2(i,j+12,k+12)=e−(σy*/μ0)ΔtHxzn−1/2(i,j+12,k+12)+1σy*Δz(1−e−(σy*/μ0)Δt)×{[Eyxn(i,j+1,k+12)−Eyxn(i,j,k+12)]⁢−[Eyzn(i,j+1,k+12)−Eyzn(i,j,k+12)]}.
The reflection of the electromagnetic wave by
the PML depended on the conductivity of each PML layer, the number of PML
layers, and the angle of incidence (with respect to the PML interface) of the incident
wave. Thus, to apply a specific PML algorithm into FDTD computations, the electric
conductivity of each layer, σi, the fictitious magnetic conductivity
of each layer, σi*, and the number of layers, NPML,
must be determined before the computations. Usually, given a required
reflection coefficient at the outer layer with normal
incidence (0° incidence angle), R(0), σi, and σi* are obtained
by [2, 3]
(10)σi=(pip)mσmax⁡,σi*=(pip)mσmax⁡*,
where p is the thickness of the PML, pi is the distance from the inner side (towards the coil) of the boundary to ith layer of PML, m is a coefficient in the range of 1 to 5, σmax⁡ and
σmax⁡* are the electric conductivity and magnetic conductivity of the
outer layer of PML respectively, which are given by
(11)σmax⁡=(m+1)ε0c02pln⁡(1R(0)),σmax*=ε0μ0σmax⁡,
where ε0 and μ0 are the permittivity and permeability in free space, respectively.According to Berenger [2, 3], an m=2 and an
R(0)=10−4 are optimal for computations using the PML algorithm. We
therefore used these values for our calculations. Generally, the reflection at the
PML interface increased with the incident angle (i.e., with deviation away from
a normal incidence). The PML was designed to absorb (with infinitesimal spatial
and temporal steps) nearly all of the electromagnetic waves having a normal incidence
(0°). However, the angle of incidence on the PML interface for RF coil
calculations is difficult to define, as radiations and reflections from the
coil or load occur at all possible angles. For practical purposes, we assumed that
the incident angle was a function of the thickness of the PML, “p,” and the distance
from the coil structure to the PML, “s” (see Figure 2). We therefore use “p”
and “s” to assess the reflection of the PML. For quantitative measurements, we
used “Np,” the layers of the PML, to denote the thickness of the PML,
“p.” Each layer was modeled by one grid. The distance “s” was also measured in
grids “Ns.”The schematic view of the computational
domain. (a) Longitudinal plane (xoz-plane in
Figure 1(b)). (b) Transverse plane
(xoy-plane in Figure 1(b)). “D” homogenous region of B1-field
(shaded); “p”: the thickness of the PML; “s”: the distance from PML to the closest
point on the coil geometry; “d”: the diameter of the coil; “l”: the length of
the coil; “dh”: the diameter of region D; “lh”: the
length of region D. The dark solid rectangle and circle denote the outline of
the coil in the longitudinal plane and transverse plane, respectively.
(a)(b)Usually, theB1 field of a
well-tuned empty birdcage coil is homogenous within a cylindrical region
located in the center of the coil. This cylindrical region, denoted by “D,”
together with the diameter “dh” and the length “lh,” is illustrated in Figure 2.


#### 2.3.2. PML Assessment Criteria

While varying the parameters of the PML, we
evaluated its cumulative error by comparing the measured
and calculatedB1 field in their resonance patterns, homogeneity, and stability. We deem a PML satisfactory
if the results of calculation met all of the following criteria.1. Resonance Pattern
For the comparison of resonance patterns, the
number of calculated modes should be equal to the number of measured modes, and
the maximum difference of the frequencies at each mode should be less than the
frequency resolution of the calculated resonance pattern (here, 0.5 MHz). That is(C1.a)Nc=Nm,(C1.b)dfimax⁡=(fci−fmi)max⁡<0.5⁢MHz,{i=0,1,…,Nm},
where 							Nc   is the number of the calculated modes;
Nm is the number of the measured modes;
dfi is the difference between the
calculated and measured frequency at ith mode; 							fci is the calculated resonance frequency at
ith mode;
fmi is the measured resonance frequency at ith mode.2.B1 Field Homogeneity
We identified a region as homogenous region
if the maximum relative deviation of theB1 field there was <10%.
This criterion, when applied to the birdcage coil, yielded
(C2)Drmax⁡=(B1c−B1¯c)max⁡B1c¯<10%,{B1c∈D},
where 				Dr⁢max⁡ is the maximum relative deviation;
B1c is the calculated B1 field;
B1¯c is the calculated mean B1 field;
D
is the cylindrical region (see Figure 2).3.B1 Field Stability
The calculatedB1 field should be
stable after implementation of a certain number of time steps. We formulated
the following criterion for this requirement:
(C3)(B1cpj−B1cpk)max⁡B1cp¯<10%,{B1c∈D;Nb<j,k<Nt},
where 				B1cpj and
B1cpk are the jth and kth peak value of the
calculated B1-field, respectively;
B1cp¯ is the mean of the peak value of the
calculated B1-field;
Nb is
the time steps at which the calculated B1-field starts to be stable;
Nt is
the total time steps of calculation.


## 2.1. Experimental Measurements

Our experimental head coil was constructed based on the high-pass
birdcage design [14] for a typical adult head. The coil had 16 rungs connected
by 2 end-rings, 32 capacitors of 13.6 pF, a diameter of 292 mm, and a length of
356 mm. The rungs and rings were made of copper strips having a width of 8 mm and
a thickness of 30 um (see Figure 1).(a) Structure of a high-pass birdcage head coil
showing a schematic view of the coil and (b) coronal and sagittal planes of the
coil, where “l”: length of the coil; “d”:
diameter of the coil; “C”: capacitor; “A”
and “B”: RF feeding points; “O”:
center of the coil.
(a)(b)We obtained the resonance pattern of the coil by measuring
the S-parameters using an Agilent 4395A network/spectrum/impedance analyzer. The
coil worked in quadrature mode [14]
with a 90° hybrid coupler
connected between the analyzer and the coil’s excitation ports. We acquired 1H
images from a phantom loaded in the coil using a GE 3 T Signa scanner and a
gradient echo pulse sequence having a 90° flip angle, 500 milliseconds
repetition time (RT]), and 10 milliseconds echo time (ET).The phantom consisted of a spherical plastic container filled
with purified water and NaCl in order to approximately simulate the average
dielectric constant (70ε0) and conductivity (0.57 S/m) of the human
brain at 128 MHz [9, 15]. The corresponding wavelength inside the phantom was
280 mm, which was 95.6% of the diameter and 78.6% of the length of the actual coil.


## 2.2. Computational Model

In simulating
the RF coil, we first used FDTD with Berenger’s PML as the boundary condition
to model the actual coil loaded with the phantom. We then computed the
resonance pattern andB1 field of the actual coil using the modeled
coil. We finally evaluated performance of the PML by comparing the calculated
results with the experimental ones. Our FDTD model and computing procedures are
described below.


### 2.2.1. FDTD Model

Our computer code for FDTD was developed in-house. The FDTD model
had 3D grids of 2 mm definition in each of thex-, y-, and z-directions. The
z-axis of the simulation was parallel to the longitudinal axis of the coil. Although
the actual coil and the phantom were modeled according to their specific
physical sizes, some of the modeled dimensions deviated from those of the real
coil because the dimensions of some of the real components were smaller than the
definition of our grid. For example, the copper strips were represented by one
voxel (2 mm) in thickness rather than by their actual thickness of 30 microns. Similarly,
the lumped capacitors were represented by one voxel and modeled by the
following equation [16]:
(1)Ixn+1/2=CΔxExi,j,kn+1−Exi,j,kn−12Δt,
where Ix was the current flowing
though the capacitor, C was the capacitance of the capacitor, and x was the
direction along which the capacitor aligned. Note that “x” was substituted with
“y” or “z” when the capacitor aligned along the respective directions.For stable time marching, the time step length,Δt must satisfy the Courant-Friedrichs-Lewy (CFL)
stability criterion [1]:
(2)Δt<cmax⁡1/Δx2+1/Δy2+1/Δz2,
where cmax⁡ is the maximum propagation velocity of the electromagnetic waves within the
simulated region. In this case, Δx=Δy=Δz=2⁢mm,
and cmax⁡, was the velocity
of light. The corresponding maximum  Δt was 3.789 picoseconds.
According to the Nyquist criterion for sampling, the associated frequency resolution
(Δf) was given as
(3)Δf=1(NΔt),
where N was the number of time steps (iterations).
For a frequency resolution of 1 MHz or less, N should be greater than 263922.
We selected N as 262144 (= 218) for convenience of calculating the fast
Fourier transformation (FFT).


### 2.2.2. Computation of Resonance Pattern

We
computed the resonance pattern of the coil by applying the FFT to its calculated
magnetic field. In so doing, we first excited the coil with an excitation signal
that contained a sufficient number of frequencies such that the calculated
electromagnetic fields of the coil contained all frequencies within its bandwidth.
Typically, the resonance frequency of the highest mode of the birdcage
coil for1H at 3 T is <500 MHz. We, therefore, chose a Gaussian
function as the excitation signal having a −3 dB frequency of 500 MHz:
(4)s(t)=e−(t/2T0−3)2,
where T0=1/ω−3⁢dB=1/(109π).We then applied this excitation signal to a capacitor of the coil and
iteratively computed the electromagnetic field of the coil. We saved the
calculatedB1 field for every time step at sampling points on the
surface and at the center of the coil (see Figure 1(b)). We then applied FFT to
the saved B1 field to obtain the frequency spectra (i.e., the resonance
pattern) of the coil at the related sampling points.A correctly calculated resonance pattern should agree with the real coil’s
actual resonance pattern provided that the frequency deviations between the
calculated and measured patterns are below the resolution of the calculated
resonance pattern and therefore are considered negligible. Otherwise, we must adjust the parameters
of the PML or the parameters of the modeled coil and recompute the resonance
pattern.


### 2.2.3. Computation of theB1 Field

After calculating the coil’s resonance
pattern, we then computed theB1 field of the coil at the Larmor
frequency of 1H at 3 T. To achieve quadrature excitation, we,
respectively, applied two sine waves centered at the Larmor frequency to feed two
capacitors 90° azimuthally apart from each other (points A and B in Figure 1(b)). The two sine wave excitations having a 90° difference in their
phases were presented as
(5a)s1(t)=(1−e−(ω0t/5)2)sin⁡(2πf0t),(5b)s2(t)=(1−e−(ω0t/5)2)sin⁡(2πf0t−π/2),
where f0 is the Larmor frequency (127.72 MHz for 1H
at 3 T). The introduction of the coefficient (1−e−(ω0t/5)2) in the excitation signal accelerated
stabilization of the calculated B1 field [10].


## 2.2.1. FDTD Model

Our computer code for FDTD was developed in-house. The FDTD model
had 3D grids of 2 mm definition in each of thex-, y-, and z-directions. The
z-axis of the simulation was parallel to the longitudinal axis of the coil. Although
the actual coil and the phantom were modeled according to their specific
physical sizes, some of the modeled dimensions deviated from those of the real
coil because the dimensions of some of the real components were smaller than the
definition of our grid. For example, the copper strips were represented by one
voxel (2 mm) in thickness rather than by their actual thickness of 30 microns. Similarly,
the lumped capacitors were represented by one voxel and modeled by the
following equation [16]:
(1)Ixn+1/2=CΔxExi,j,kn+1−Exi,j,kn−12Δt,
where Ix was the current flowing
though the capacitor, C was the capacitance of the capacitor, and x was the
direction along which the capacitor aligned. Note that “x” was substituted with
“y” or “z” when the capacitor aligned along the respective directions.For stable time marching, the time step length,Δt must satisfy the Courant-Friedrichs-Lewy (CFL)
stability criterion [1]:
(2)Δt<cmax⁡1/Δx2+1/Δy2+1/Δz2,
where cmax⁡ is the maximum propagation velocity of the electromagnetic waves within the
simulated region. In this case, Δx=Δy=Δz=2⁢mm,
and cmax⁡, was the velocity
of light. The corresponding maximum  Δt was 3.789 picoseconds.
According to the Nyquist criterion for sampling, the associated frequency resolution
(Δf) was given as
(3)Δf=1(NΔt),
where N was the number of time steps (iterations).
For a frequency resolution of 1 MHz or less, N should be greater than 263922.
We selected N as 262144 (= 218) for convenience of calculating the fast
Fourier transformation (FFT).


## 2.2.2. Computation of Resonance Pattern

We
computed the resonance pattern of the coil by applying the FFT to its calculated
magnetic field. In so doing, we first excited the coil with an excitation signal
that contained a sufficient number of frequencies such that the calculated
electromagnetic fields of the coil contained all frequencies within its bandwidth.
Typically, the resonance frequency of the highest mode of the birdcage
coil for1H at 3 T is <500 MHz. We, therefore, chose a Gaussian
function as the excitation signal having a −3 dB frequency of 500 MHz:
(4)s(t)=e−(t/2T0−3)2,
where T0=1/ω−3⁢dB=1/(109π).We then applied this excitation signal to a capacitor of the coil and
iteratively computed the electromagnetic field of the coil. We saved the
calculatedB1 field for every time step at sampling points on the
surface and at the center of the coil (see Figure 1(b)). We then applied FFT to
the saved B1 field to obtain the frequency spectra (i.e., the resonance
pattern) of the coil at the related sampling points.A correctly calculated resonance pattern should agree with the real coil’s
actual resonance pattern provided that the frequency deviations between the
calculated and measured patterns are below the resolution of the calculated
resonance pattern and therefore are considered negligible. Otherwise, we must adjust the parameters
of the PML or the parameters of the modeled coil and recompute the resonance
pattern.


## 2.2.3. Computation of theB1 Field

After calculating the coil’s resonance
pattern, we then computed theB1 field of the coil at the Larmor
frequency of 1H at 3 T. To achieve quadrature excitation, we,
respectively, applied two sine waves centered at the Larmor frequency to feed two
capacitors 90° azimuthally apart from each other (points A and B in Figure 1(b)). The two sine wave excitations having a 90° difference in their
phases were presented as
(5a)s1(t)=(1−e−(ω0t/5)2)sin⁡(2πf0t),(5b)s2(t)=(1−e−(ω0t/5)2)sin⁡(2πf0t−π/2),
where f0 is the Larmor frequency (127.72 MHz for 1H
at 3 T). The introduction of the coefficient (1−e−(ω0t/5)2) in the excitation signal accelerated
stabilization of the calculated B1 field [10].


## 2.3. PML

### 2.3.1. PML Parameters

When implementing Berenger’s PML [2, 3], layers
of artificial perfect electric conductor (PEC) were introduced at the boundaries
of the computational domain to absorb the waves traveling outwards (i.e., to satisfy
the radiation condition). The conductivity of each layer increased gradually
from inside (near the coil’s geometry) to outside (away from the coil’s
geometry). Thus, the outgoing waves attenuated to negligible levels as they
traveled through the PML, then reflected at the PEC walls located at the outer
boundaries of the computational domain, and then traveled through the PML back
to the geometry of the coil.To compute the fields in the presence of Berenger’s
PML, each component of the electromagnetic fields must be divided into two
subcomponents. For example, the electric field componentEx was
divided into Exy and Exz:
(6)Ex=Exy+Exz,
where
(7a)Exyn+1(i+12,j,k)=e−(σy/ε0)ΔtExyn(i+12,j,k)+1σyΔy(1−e−(σy/ε0)Δt)×{[Hzxn+1/2(i+12,j+12,k)−Hzxn+1/2(i+12,j−12,k)]⁢−[Hzyn+1/2(i+12,j,k+12)−Hzyn+1/2(i+12,j,k−12)]},(7b)Exzn+1(i+12,j,k)=e−(σz/ε0)ΔtExzn(i+12,j,k)+1σzΔz(1−e−(σz/ε0)Δt)×{[Hyxn+1/2(i+12,j+12,k)−Hyxn+1/2(i+12,j−12,k)]⁢−[Hyzn+1/2(i+12,j,k+12)−Hyzn+1/2(i+12,j,k−12)]}.
Similarly, the magnetic field component Hx was divided
into Hxy and Hxz:
(8)Hx=Hxy+Hxz,
where
(9a)Hxyn+1/2(i,j+12,k+12)=e−(σy*/μ0)ΔtHxyn−1/2(i,j+12,k+12)+1σy*Δy(1−e−(σy*/μ0)Δt)×{[Ezxn(i,j+1,k+12)−Ezxn(i,j,k+12)]⁢−[Ezyn(i,j+1,k+12)−Ezyn(i,j,k+12)]},(9b)Hxzn+1/2(i,j+12,k+12)=e−(σy*/μ0)ΔtHxzn−1/2(i,j+12,k+12)+1σy*Δz(1−e−(σy*/μ0)Δt)×{[Eyxn(i,j+1,k+12)−Eyxn(i,j,k+12)]⁢−[Eyzn(i,j+1,k+12)−Eyzn(i,j,k+12)]}.
The reflection of the electromagnetic wave by
the PML depended on the conductivity of each PML layer, the number of PML
layers, and the angle of incidence (with respect to the PML interface) of the incident
wave. Thus, to apply a specific PML algorithm into FDTD computations, the electric
conductivity of each layer, σi, the fictitious magnetic conductivity
of each layer, σi*, and the number of layers, NPML,
must be determined before the computations. Usually, given a required
reflection coefficient at the outer layer with normal
incidence (0° incidence angle), R(0), σi, and σi* are obtained
by [2, 3]
(10)σi=(pip)mσmax⁡,σi*=(pip)mσmax⁡*,
where p is the thickness of the PML, pi is the distance from the inner side (towards the coil) of the boundary to ith layer of PML, m is a coefficient in the range of 1 to 5, σmax⁡ and
σmax⁡* are the electric conductivity and magnetic conductivity of the
outer layer of PML respectively, which are given by
(11)σmax⁡=(m+1)ε0c02pln⁡(1R(0)),σmax*=ε0μ0σmax⁡,
where ε0 and μ0 are the permittivity and permeability in free space, respectively.According to Berenger [2, 3], an m=2 and an
R(0)=10−4 are optimal for computations using the PML algorithm. We
therefore used these values for our calculations. Generally, the reflection at the
PML interface increased with the incident angle (i.e., with deviation away from
a normal incidence). The PML was designed to absorb (with infinitesimal spatial
and temporal steps) nearly all of the electromagnetic waves having a normal incidence
(0°). However, the angle of incidence on the PML interface for RF coil
calculations is difficult to define, as radiations and reflections from the
coil or load occur at all possible angles. For practical purposes, we assumed that
the incident angle was a function of the thickness of the PML, “p,” and the distance
from the coil structure to the PML, “s” (see Figure 2). We therefore use “p”
and “s” to assess the reflection of the PML. For quantitative measurements, we
used “Np,” the layers of the PML, to denote the thickness of the PML,
“p.” Each layer was modeled by one grid. The distance “s” was also measured in
grids “Ns.”The schematic view of the computational
domain. (a) Longitudinal plane (xoz-plane in
Figure 1(b)). (b) Transverse plane
(xoy-plane in Figure 1(b)). “D” homogenous region of B1-field
(shaded); “p”: the thickness of the PML; “s”: the distance from PML to the closest
point on the coil geometry; “d”: the diameter of the coil; “l”: the length of
the coil; “dh”: the diameter of region D; “lh”: the
length of region D. The dark solid rectangle and circle denote the outline of
the coil in the longitudinal plane and transverse plane, respectively.
(a)(b)Usually, theB1 field of a
well-tuned empty birdcage coil is homogenous within a cylindrical region
located in the center of the coil. This cylindrical region, denoted by “D,”
together with the diameter “dh” and the length “lh,” is illustrated in Figure 2.


### 2.3.2. PML Assessment Criteria

While varying the parameters of the PML, we
evaluated its cumulative error by comparing the measured
and calculatedB1 field in their resonance patterns, homogeneity, and stability. We deem a PML satisfactory
if the results of calculation met all of the following criteria.1. Resonance Pattern
For the comparison of resonance patterns, the
number of calculated modes should be equal to the number of measured modes, and
the maximum difference of the frequencies at each mode should be less than the
frequency resolution of the calculated resonance pattern (here, 0.5 MHz). That is(C1.a)Nc=Nm,(C1.b)dfimax⁡=(fci−fmi)max⁡<0.5⁢MHz,{i=0,1,…,Nm},
where 							Nc   is the number of the calculated modes;
Nm is the number of the measured modes;
dfi is the difference between the
calculated and measured frequency at ith mode; 							fci is the calculated resonance frequency at
ith mode;
fmi is the measured resonance frequency at ith mode.2.B1 Field Homogeneity
We identified a region as homogenous region
if the maximum relative deviation of theB1 field there was <10%.
This criterion, when applied to the birdcage coil, yielded
(C2)Drmax⁡=(B1c−B1¯c)max⁡B1c¯<10%,{B1c∈D},
where 				Dr⁢max⁡ is the maximum relative deviation;
B1c is the calculated B1 field;
B1¯c is the calculated mean B1 field;
D
is the cylindrical region (see Figure 2).3.B1 Field Stability
The calculatedB1 field should be
stable after implementation of a certain number of time steps. We formulated
the following criterion for this requirement:
(C3)(B1cpj−B1cpk)max⁡B1cp¯<10%,{B1c∈D;Nb<j,k<Nt},
where 				B1cpj and
B1cpk are the jth and kth peak value of the
calculated B1-field, respectively;
B1cp¯ is the mean of the peak value of the
calculated B1-field;
Nb is
the time steps at which the calculated B1-field starts to be stable;
Nt is
the total time steps of calculation.


## 2.3.1. PML Parameters

When implementing Berenger’s PML [2, 3], layers
of artificial perfect electric conductor (PEC) were introduced at the boundaries
of the computational domain to absorb the waves traveling outwards (i.e., to satisfy
the radiation condition). The conductivity of each layer increased gradually
from inside (near the coil’s geometry) to outside (away from the coil’s
geometry). Thus, the outgoing waves attenuated to negligible levels as they
traveled through the PML, then reflected at the PEC walls located at the outer
boundaries of the computational domain, and then traveled through the PML back
to the geometry of the coil.To compute the fields in the presence of Berenger’s
PML, each component of the electromagnetic fields must be divided into two
subcomponents. For example, the electric field componentEx was
divided into Exy and Exz:
(6)Ex=Exy+Exz,
where
(7a)Exyn+1(i+12,j,k)=e−(σy/ε0)ΔtExyn(i+12,j,k)+1σyΔy(1−e−(σy/ε0)Δt)×{[Hzxn+1/2(i+12,j+12,k)−Hzxn+1/2(i+12,j−12,k)]⁢−[Hzyn+1/2(i+12,j,k+12)−Hzyn+1/2(i+12,j,k−12)]},(7b)Exzn+1(i+12,j,k)=e−(σz/ε0)ΔtExzn(i+12,j,k)+1σzΔz(1−e−(σz/ε0)Δt)×{[Hyxn+1/2(i+12,j+12,k)−Hyxn+1/2(i+12,j−12,k)]⁢−[Hyzn+1/2(i+12,j,k+12)−Hyzn+1/2(i+12,j,k−12)]}.
Similarly, the magnetic field component Hx was divided
into Hxy and Hxz:
(8)Hx=Hxy+Hxz,
where
(9a)Hxyn+1/2(i,j+12,k+12)=e−(σy*/μ0)ΔtHxyn−1/2(i,j+12,k+12)+1σy*Δy(1−e−(σy*/μ0)Δt)×{[Ezxn(i,j+1,k+12)−Ezxn(i,j,k+12)]⁢−[Ezyn(i,j+1,k+12)−Ezyn(i,j,k+12)]},(9b)Hxzn+1/2(i,j+12,k+12)=e−(σy*/μ0)ΔtHxzn−1/2(i,j+12,k+12)+1σy*Δz(1−e−(σy*/μ0)Δt)×{[Eyxn(i,j+1,k+12)−Eyxn(i,j,k+12)]⁢−[Eyzn(i,j+1,k+12)−Eyzn(i,j,k+12)]}.
The reflection of the electromagnetic wave by
the PML depended on the conductivity of each PML layer, the number of PML
layers, and the angle of incidence (with respect to the PML interface) of the incident
wave. Thus, to apply a specific PML algorithm into FDTD computations, the electric
conductivity of each layer, σi, the fictitious magnetic conductivity
of each layer, σi*, and the number of layers, NPML,
must be determined before the computations. Usually, given a required
reflection coefficient at the outer layer with normal
incidence (0° incidence angle), R(0), σi, and σi* are obtained
by [2, 3]
(10)σi=(pip)mσmax⁡,σi*=(pip)mσmax⁡*,
where p is the thickness of the PML, pi is the distance from the inner side (towards the coil) of the boundary to ith layer of PML, m is a coefficient in the range of 1 to 5, σmax⁡ and
σmax⁡* are the electric conductivity and magnetic conductivity of the
outer layer of PML respectively, which are given by
(11)σmax⁡=(m+1)ε0c02pln⁡(1R(0)),σmax*=ε0μ0σmax⁡,
where ε0 and μ0 are the permittivity and permeability in free space, respectively.According to Berenger [2, 3], an m=2 and an
R(0)=10−4 are optimal for computations using the PML algorithm. We
therefore used these values for our calculations. Generally, the reflection at the
PML interface increased with the incident angle (i.e., with deviation away from
a normal incidence). The PML was designed to absorb (with infinitesimal spatial
and temporal steps) nearly all of the electromagnetic waves having a normal incidence
(0°). However, the angle of incidence on the PML interface for RF coil
calculations is difficult to define, as radiations and reflections from the
coil or load occur at all possible angles. For practical purposes, we assumed that
the incident angle was a function of the thickness of the PML, “p,” and the distance
from the coil structure to the PML, “s” (see Figure 2). We therefore use “p”
and “s” to assess the reflection of the PML. For quantitative measurements, we
used “Np,” the layers of the PML, to denote the thickness of the PML,
“p.” Each layer was modeled by one grid. The distance “s” was also measured in
grids “Ns.”The schematic view of the computational
domain. (a) Longitudinal plane (xoz-plane in
Figure 1(b)). (b) Transverse plane
(xoy-plane in Figure 1(b)). “D” homogenous region of B1-field
(shaded); “p”: the thickness of the PML; “s”: the distance from PML to the closest
point on the coil geometry; “d”: the diameter of the coil; “l”: the length of
the coil; “dh”: the diameter of region D; “lh”: the
length of region D. The dark solid rectangle and circle denote the outline of
the coil in the longitudinal plane and transverse plane, respectively.
(a)(b)Usually, theB1 field of a
well-tuned empty birdcage coil is homogenous within a cylindrical region
located in the center of the coil. This cylindrical region, denoted by “D,”
together with the diameter “dh” and the length “lh,” is illustrated in Figure 2.


## 2.3.2. PML Assessment Criteria

While varying the parameters of the PML, we
evaluated its cumulative error by comparing the measured
and calculatedB1 field in their resonance patterns, homogeneity, and stability. We deem a PML satisfactory
if the results of calculation met all of the following criteria.1. Resonance Pattern
For the comparison of resonance patterns, the
number of calculated modes should be equal to the number of measured modes, and
the maximum difference of the frequencies at each mode should be less than the
frequency resolution of the calculated resonance pattern (here, 0.5 MHz). That is(C1.a)Nc=Nm,(C1.b)dfimax⁡=(fci−fmi)max⁡<0.5⁢MHz,{i=0,1,…,Nm},
where 							Nc   is the number of the calculated modes;
Nm is the number of the measured modes;
dfi is the difference between the
calculated and measured frequency at ith mode; 							fci is the calculated resonance frequency at
ith mode;
fmi is the measured resonance frequency at ith mode.2.B1 Field Homogeneity
We identified a region as homogenous region
if the maximum relative deviation of theB1 field there was <10%.
This criterion, when applied to the birdcage coil, yielded
(C2)Drmax⁡=(B1c−B1¯c)max⁡B1c¯<10%,{B1c∈D},
where 				Dr⁢max⁡ is the maximum relative deviation;
B1c is the calculated B1 field;
B1¯c is the calculated mean B1 field;
D
is the cylindrical region (see Figure 2).3.B1 Field Stability
The calculatedB1 field should be
stable after implementation of a certain number of time steps. We formulated
the following criterion for this requirement:
(C3)(B1cpj−B1cpk)max⁡B1cp¯<10%,{B1c∈D;Nb<j,k<Nt},
where 				B1cpj and
B1cpk are the jth and kth peak value of the
calculated B1-field, respectively;
B1cp¯ is the mean of the peak value of the
calculated B1-field;
Nb is
the time steps at which the calculated B1-field starts to be stable;
Nt is
the total time steps of calculation.


## 3. Results and Discussion

We found that both the resonance pattern and homogeneity
of the calculatedB1-field agreed with the measurements of the real
coil when using an optimal PML but not when using nonoptimal PMLs.


### 3.1. Resonance Pattern

The
measured resonance pattern had eight distinct resonance modes (see Figure3(a)).
Mode 0 had the highest frequency at 180.2 MHz. Mode 1 (the operational mode, marked
as H1) had a frequency of 127.72 MHz, which is the Larmor frequency
of 1H at 3 T. The other modes had frequencies ranging between 57.8
and 96.4 MHz.(a) MeasuredS11 of the coil
(obtained with a network analyzer). H1 is the mode of interest (mode 1); (b) calculated
resonance pattern with Np=5 and Ns=5 grids. Only one mode
resonated at 174 MHz; (c) calculated resonance pattern with Np=10 and
Ns=8 grids. The number of modes and resonance frequencies is  in excellent agreement with the measured S11;
(d) calculated resonance pattern with Np=10 and Ns=1
grids. Mode 0 and mode 1 are split (e.g., mode 1 split into two peaks: H1 and
H1′.)
(a)(b)(c)(d)Intuitively, a good PML should be sufficiently
thick to attenuate outgoing waves without contaminating the actual
electromagnetic fields within the coil. A good PML also should be sufficiently
distant from the coil’s geometry to yield the smallest possible incidence angles
of the outgoing waves at the PML interface. By varying the thickness of the PML,Np and the distance from the PML to the closest surface of the coil, Ns, we obtained the calculated
resonance pattern in which frequencies agreed with the measured resonance
pattern (see Figure 3(c)). Because the resolution of the calculated resonance
pattern was 0.5 MHz after 262144 time steps, the frequency of each calculated
and measured modes inevitably deviated. However, the maximum deviation between
the calculated and measured frequencies did not exceed the resolution of the calculated
frequency, indicating that this model provided highly accurate estimates of the
measured frequencies.We
found that to obtain the correct resonance pattern for this particular coil,
the number of layers of the PML,Np, must be greater than 8 and the number
of grids between the closest point on the coil’s geometry and the PML interface,
Ns, must be greater than 5. Any PML having an Np<8 or an Ns<5 produced a resonance at a single frequency (see Figure 3(b)) or a split
of the resonance modes (see Figure 3(d)). Based on these results for the computation
of the coil’s resonance pattern, we advise use of a PML having a thickness of
at least 10 layers and a distance of at least 8 grids between the closest point
on the coil’s geometry and the PML interface. Note that these observations
obviously are valid for the particular spatial and temporal steps that we
examined. Nonetheless, splitting of the resonance modes of an RF coil in
simulations but not in the bench testing could be caused by a poor choice of
PML parameters.


### 3.2.B1 Homogeneity and Stability

To evaluateB1 homogeneity and stability,
we first computed the B1 field of the coil using a PML with higher Np and Ns values and then we repeated the computations using a PML with
lower Np and Ns values. Comparing the two computed
results demonstrated that when the PML with higher Np
and Ns values increased the amplitude of the calculated B1-field gradually during
the first 40 000 time steps (corresponding to about 40 cycles of
the B1 field at 127.72 MHz) and then stabilized thereafter (see Figure 4(a)). Using the
PML with lower Np and Ns values, however, increased the
amplitude of the calculated B1-field monotonically and without bound
(see Figure 4(b)).CalculatedB1-field at the center of the coil. (a)
PML with Np=10 and Ns=10. The B1-field is
stable after 40 000 time steps; (b) PML with Np=5 and Ns=5.
The B1-field increases monotonically even after 100 000 time steps.
(a)(b)These findings indicate that we obtained an optimal
PML by iteratively calculating and evaluating theB1 field of the
coil for either 40 or 120 cycles and by simultaneously varying the thickness of
the PML and the distance between the PML interface and the closest point on the
coil’s geometry. We let Np equal 5, 8, 10, or 12 while increasing Ns from 1 to 20 for each value of Np. While performing these numerical
studies of the effects of varying Np and Ns, we
calculated the coil’s B1 field and evaluated its homogeneity
according to our previously stated criterion C2 and C3 for each PML setting. These
assessments revealed that the homogeneity of the calculated B1 field
improved with increases of both Np and Ns (see Figure 5).
We noted that the calculated deviations were especially sensitive to values of
Np and Ns that were below 8 and 5, respectively, for 40-cycle
computations, whereas increasing Np and Ns beyond 8 and
6, respectively, did not noticeably improve homogeneity of the calculated B1 field (see Figure 5(a)). When using 120 cycles, the minimum optimal Np and Ns increased to 10 layers and 12 grids (see Figure 5(b)).Maximum deviations of the calculatedB1-field within
a circular region positioned in axial plane at the center of the coil. This
region has a diameter that is 80% of the coil’s diameter. The results are
presented at 127.72 MHz, the (a) Larmor frequency of 1H at 3 T after
iterations of 40 cycles and (b) after iterations of 120 cycles.
(a)(b)Because computational time is proportional to
the number of grids, which in turn is directly related to the sum ofNp and Ns to the third power, Np and Ns should be
minimized to reduce the computational time, but only so far as to retain the
required degree of computational accuracy. Thus for the calculation of a 40-cycle
B1 field, an optimal PML should be 8 layers thick and at a distance
of 5 grids from the closest point on the coil’s geometry (see Figure 5(a)). For
the calculation of a 120-cycle B1 field, an optimal PML should be at least 10 layers thick and at a distance of 12
grids from the coil’s geometry (see Figure 5(b)).Additionally, we found that a well-optimized
PML yielded a homogeneous calculatedB1 field (see Figure 6(b)) that
agreed well with the measured field (see Figure 6(a)). A poorly optimized PML,
however, yielded an inhomogeneous calculated B1 field both inside
and outside the coil (see Figure 6(c))). Furthermore, a PML that satisfied the
homogeneity criterion after 40 000 time steps can deteriorate when using a
larger number of time steps (see Figure 6(d)), because the reflections on the PML
can accumulate to unacceptable levels with the increasing number of iterations.Measurement of homogeneity ofB1 field. (a) Acquired MR image
of an oil phantom loaded within the coil. (b) Calculated B1-map
after 40-cycle iterations with Np=10 and Ns=10. The
homogeneity of B1-map is good inside, on the surface of, and outside
of the coil. (c) Calculated B1-map after 40-cycle iterations with Np=10
and Ns=1. The homogeneity of B1-map is unacceptable. (d)
Calculated B1-map after 120-cycle iterations with Np=10
and Ns=10. The intensities of the B1 field at the
left-bottom corner and right-top corner are higher than those at other corners.
(a)(b)(c)(d)


### 3.3. Simulation with a Heterogeneous Phantom

The above findings were based on tests using a
homogenous phantom. A human head, however, consists of various tissues, each of
which differs in their dielectric and conductivity properties (see Table1)
[10, 13, 15]. To assess the influences of these heterogeneous tissues on
the accuracy of the PML, we simulated a realistic human head using a
heterogeneous phantom and then computed the resonance pattern and B1 fields of the coil in which it was placed. The outer layer of grids of the heterogeneous
phantom was modeled as skin (2 mm thick), three layers of grids inside the skin
layer were modeled as skull (6 mm thick), and another three layers next to
skull were modeled as cerebrospinal fluid (CSF). The remaining center portion of
the phantom was modeled as brain.Table 1
The dielectric and
conductivity of tissues inside a human head at 128 MHz.
Tissue typeεrσ (s/m)Skin630.86Skull250.05CSF1030.7Brain700.57The calculated results showed that when using
the same PML parameters, (a) the resonance frequency at each mode of the
calculated resonance pattern kept unchanged, even though the calculatedQ-factor of the coil loaded with the heterogeneous phantom compared with that
of the homogeneous phantom decreased by 4.6% at the operating frequency of
128 MHz; (b) the deviation of the calculated B1 field of the coil loaded
with the heterogeneous phantom compared with that loaded with homogeneous
phantom increased by 1.4%. These changes in both the Q-factor and deviation of the
B1 field were caused by differences in heterogeneity of the phantom rather
than by the PML used in the computations. Thus, we believe that the parameters
of a PML optimized for homogeneous phantoms also apply to heterogeneous
phantoms.


## 3.1. Resonance Pattern

The
measured resonance pattern had eight distinct resonance modes (see Figure3(a)).
Mode 0 had the highest frequency at 180.2 MHz. Mode 1 (the operational mode, marked
as H1) had a frequency of 127.72 MHz, which is the Larmor frequency
of 1H at 3 T. The other modes had frequencies ranging between 57.8
and 96.4 MHz.(a) MeasuredS11 of the coil
(obtained with a network analyzer). H1 is the mode of interest (mode 1); (b) calculated
resonance pattern with Np=5 and Ns=5 grids. Only one mode
resonated at 174 MHz; (c) calculated resonance pattern with Np=10 and
Ns=8 grids. The number of modes and resonance frequencies is  in excellent agreement with the measured S11;
(d) calculated resonance pattern with Np=10 and Ns=1
grids. Mode 0 and mode 1 are split (e.g., mode 1 split into two peaks: H1 and
H1′.)
(a)(b)(c)(d)Intuitively, a good PML should be sufficiently
thick to attenuate outgoing waves without contaminating the actual
electromagnetic fields within the coil. A good PML also should be sufficiently
distant from the coil’s geometry to yield the smallest possible incidence angles
of the outgoing waves at the PML interface. By varying the thickness of the PML,Np and the distance from the PML to the closest surface of the coil, Ns, we obtained the calculated
resonance pattern in which frequencies agreed with the measured resonance
pattern (see Figure 3(c)). Because the resolution of the calculated resonance
pattern was 0.5 MHz after 262144 time steps, the frequency of each calculated
and measured modes inevitably deviated. However, the maximum deviation between
the calculated and measured frequencies did not exceed the resolution of the calculated
frequency, indicating that this model provided highly accurate estimates of the
measured frequencies.We
found that to obtain the correct resonance pattern for this particular coil,
the number of layers of the PML,Np, must be greater than 8 and the number
of grids between the closest point on the coil’s geometry and the PML interface,
Ns, must be greater than 5. Any PML having an Np<8 or an Ns<5 produced a resonance at a single frequency (see Figure 3(b)) or a split
of the resonance modes (see Figure 3(d)). Based on these results for the computation
of the coil’s resonance pattern, we advise use of a PML having a thickness of
at least 10 layers and a distance of at least 8 grids between the closest point
on the coil’s geometry and the PML interface. Note that these observations
obviously are valid for the particular spatial and temporal steps that we
examined. Nonetheless, splitting of the resonance modes of an RF coil in
simulations but not in the bench testing could be caused by a poor choice of
PML parameters.


## 3.2.B1 Homogeneity and Stability

To evaluateB1 homogeneity and stability,
we first computed the B1 field of the coil using a PML with higher Np and Ns values and then we repeated the computations using a PML with
lower Np and Ns values. Comparing the two computed
results demonstrated that when the PML with higher Np
and Ns values increased the amplitude of the calculated B1-field gradually during
the first 40 000 time steps (corresponding to about 40 cycles of
the B1 field at 127.72 MHz) and then stabilized thereafter (see Figure 4(a)). Using the
PML with lower Np and Ns values, however, increased the
amplitude of the calculated B1-field monotonically and without bound
(see Figure 4(b)).CalculatedB1-field at the center of the coil. (a)
PML with Np=10 and Ns=10. The B1-field is
stable after 40 000 time steps; (b) PML with Np=5 and Ns=5.
The B1-field increases monotonically even after 100 000 time steps.
(a)(b)These findings indicate that we obtained an optimal
PML by iteratively calculating and evaluating theB1 field of the
coil for either 40 or 120 cycles and by simultaneously varying the thickness of
the PML and the distance between the PML interface and the closest point on the
coil’s geometry. We let Np equal 5, 8, 10, or 12 while increasing Ns from 1 to 20 for each value of Np. While performing these numerical
studies of the effects of varying Np and Ns, we
calculated the coil’s B1 field and evaluated its homogeneity
according to our previously stated criterion C2 and C3 for each PML setting. These
assessments revealed that the homogeneity of the calculated B1 field
improved with increases of both Np and Ns (see Figure 5).
We noted that the calculated deviations were especially sensitive to values of
Np and Ns that were below 8 and 5, respectively, for 40-cycle
computations, whereas increasing Np and Ns beyond 8 and
6, respectively, did not noticeably improve homogeneity of the calculated B1 field (see Figure 5(a)). When using 120 cycles, the minimum optimal Np and Ns increased to 10 layers and 12 grids (see Figure 5(b)).Maximum deviations of the calculatedB1-field within
a circular region positioned in axial plane at the center of the coil. This
region has a diameter that is 80% of the coil’s diameter. The results are
presented at 127.72 MHz, the (a) Larmor frequency of 1H at 3 T after
iterations of 40 cycles and (b) after iterations of 120 cycles.
(a)(b)Because computational time is proportional to
the number of grids, which in turn is directly related to the sum ofNp and Ns to the third power, Np and Ns should be
minimized to reduce the computational time, but only so far as to retain the
required degree of computational accuracy. Thus for the calculation of a 40-cycle
B1 field, an optimal PML should be 8 layers thick and at a distance
of 5 grids from the closest point on the coil’s geometry (see Figure 5(a)). For
the calculation of a 120-cycle B1 field, an optimal PML should be at least 10 layers thick and at a distance of 12
grids from the coil’s geometry (see Figure 5(b)).Additionally, we found that a well-optimized
PML yielded a homogeneous calculatedB1 field (see Figure 6(b)) that
agreed well with the measured field (see Figure 6(a)). A poorly optimized PML,
however, yielded an inhomogeneous calculated B1 field both inside
and outside the coil (see Figure 6(c))). Furthermore, a PML that satisfied the
homogeneity criterion after 40 000 time steps can deteriorate when using a
larger number of time steps (see Figure 6(d)), because the reflections on the PML
can accumulate to unacceptable levels with the increasing number of iterations.Measurement of homogeneity ofB1 field. (a) Acquired MR image
of an oil phantom loaded within the coil. (b) Calculated B1-map
after 40-cycle iterations with Np=10 and Ns=10. The
homogeneity of B1-map is good inside, on the surface of, and outside
of the coil. (c) Calculated B1-map after 40-cycle iterations with Np=10
and Ns=1. The homogeneity of B1-map is unacceptable. (d)
Calculated B1-map after 120-cycle iterations with Np=10
and Ns=10. The intensities of the B1 field at the
left-bottom corner and right-top corner are higher than those at other corners.
(a)(b)(c)(d)


## 3.3. Simulation with a Heterogeneous Phantom

The above findings were based on tests using a
homogenous phantom. A human head, however, consists of various tissues, each of
which differs in their dielectric and conductivity properties (see Table1)
[10, 13, 15]. To assess the influences of these heterogeneous tissues on
the accuracy of the PML, we simulated a realistic human head using a
heterogeneous phantom and then computed the resonance pattern and B1 fields of the coil in which it was placed. The outer layer of grids of the heterogeneous
phantom was modeled as skin (2 mm thick), three layers of grids inside the skin
layer were modeled as skull (6 mm thick), and another three layers next to
skull were modeled as cerebrospinal fluid (CSF). The remaining center portion of
the phantom was modeled as brain.Table 1
The dielectric and
conductivity of tissues inside a human head at 128 MHz.
Tissue typeεrσ (s/m)Skin630.86Skull250.05CSF1030.7Brain700.57The calculated results showed that when using
the same PML parameters, (a) the resonance frequency at each mode of the
calculated resonance pattern kept unchanged, even though the calculatedQ-factor of the coil loaded with the heterogeneous phantom compared with that
of the homogeneous phantom decreased by 4.6% at the operating frequency of
128 MHz; (b) the deviation of the calculated B1 field of the coil loaded
with the heterogeneous phantom compared with that loaded with homogeneous
phantom increased by 1.4%. These changes in both the Q-factor and deviation of the
B1 field were caused by differences in heterogeneity of the phantom rather
than by the PML used in the computations. Thus, we believe that the parameters
of a PML optimized for homogeneous phantoms also apply to heterogeneous
phantoms.


## 4. Conclusion

We have explored the role of Berenger’s PML as
an absorbing boundary condition in the computational characterization of RF
coils for MRI at 3 T. We presented a method that evaluates the accuracy and
efficiency of the PML for computing the resonance patterns andB1 fields
of RF coils using FDTD. We defined three criteria to evaluate the PML: (1) the deviation
of the frequencies of the resonance mode for the calculated and the measured
data, (2) the maximum relative deviation of the calculated B1 field
inside the coil, and (3) the numerical stability of the calculated B1 field
of the coil. The results demonstrated that the accuracy and efficiency of the PML
as an absorbing boundary condition closely depend on both the thickness of the
PML and the space between the PML and the coil surface. For computation of the
coil’s resonance patterns, a poor choice of PML parameters can produce either a
split in the resonance modes or a resonance that occurs only at a single mode.
For the computation of B1 field, a poor choice of PML parameters can
produce a large fluctuation of the B1 field. Additionally, PML
parameters that yield good computational accuracy for a specified number of time
steps can also yield inaccurate results with an excessive number of time steps,
implying the excess accumulation of reflections of the RF wave on the PML. Taking
into account both accuracy and efficiency for specified spatial and temporal
steps, we found that an optimal PML for both the computation of the resonance
pattern and the 40-cycle B1 field of a 3 T head coil loaded with a
phantom has a thickness of at least 8 layers and a distance of 5 grids from the
PML to the closest point on the coil’s geometry. Furthermore, the thickness will
increase to at least 10 layers and the distance to at least 12 grids for computations
of more than 40 cycles of the B1 fields.



---

*Source: 563196-2008-09-03.xml*