International Journal of Antennas and Propagation
Volume 2008 (2008), Article ID 563196, 10 pages
doi:10.1155/2008/563196
Research Article

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

Yunsuo Duan,1 Tamer S. Ibrahim,2,3 Bradley S. Peterson,1 Feng Liu,1 and Alayar Kangarlu1

1Columbia University and New York State Psychiatric Institute, New York, NY 10032, USA
2University of Pittsburgh, Pittsburgh, PA 15260, USA
3University of Oklahoma, Norman, OK 73019-0390, USA

Received 16 October 2007; Revised 9 April 2008; Accepted 5 June 2008

Academic Editor: Stuart Crozier

Copyright © 2008 Yunsuo Duan et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

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 ( 1 H ) phantom MRI images. We then computed the resonance pattern and B 1 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 1 H at 3 Tesla (T).

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 ( 𝐵 0 ) 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 𝐵 0 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 of 1H 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 𝐵 1 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 the 𝐵 1 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 𝐵 1 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).

fig1
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.

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 and 𝐵 1 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 the x-, 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 / 2 𝐸 = 𝐶 Δ 𝑥 𝑥 𝑛 + 1 𝑖 , 𝑗 , 𝑘 − 𝐸 𝑥 𝑛 − 1 𝑖 , 𝑗 , 𝑘 2 Δ 𝑡 , ( 1 ) where 𝐼 𝑥 was the current flowing though the capacitor, 𝐶 was the capacitance of the capacitor, and 𝑥 was the direction along which the capacitor aligned. Note that “ 𝑥 ” was substituted with “ 𝑦 ” or “ 𝑧 ” when the capacitor aligned along the respective directions.

For stable time marching, the time step length, Δ 𝑡 must satisfy the Courant-Friedrichs-Lewy (CFL) stability criterion [1]: 𝑐 Δ 𝑡 < m a x √ 1 / Δ 𝑥 2 + 1 / Δ 𝑦 2 + 1 / Δ 𝑧 2 , ( 2 ) where 𝑐 m a x is the maximum propagation velocity of the electromagnetic waves within the simulated region. In this case, Δ 𝑥 = Δ 𝑦 = Δ 𝑧 = 2 m m , and 𝑐 m a x , was the velocity of light. The corresponding maximum Δ 𝑡 was 3.789 picoseconds. According to the Nyquist criterion for sampling, the associated frequency resolution ( Δ 𝑓 ) was given as 1 Δ 𝑓 = ( 𝑁 Δ 𝑡 ) , ( 3 ) where 𝑁 was the number of time steps (iterations). For a frequency resolution of 1 MHz or less, 𝑁 should be greater than 263922. We selected 𝑁 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 for 1H at 3 T is <500 MHz. We, therefore, chose a Gaussian function as the excitation signal having a −3 dB frequency of 500 MHz: 𝑠 ( 𝑡 ) = 𝑒 − ( 𝑡 / 2 𝑇 0 − 3 ) 2 , ( 4 ) where 𝑇 0 = 1 / 𝜔 − 3 d B = 1 / ( 1 0 9 𝜋 ) .

We then applied this excitation signal to a capacitor of the coil and iteratively computed the electromagnetic field of the coil. We saved the calculated 𝐵 1 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 𝐵 1 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 the 𝐵 1 Field

After calculating the coil’s resonance pattern, we then computed the 𝐵 1 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 𝐴 and 𝐵 in Figure 1(b)). The two sine wave excitations having a 90° difference in their phases were presented as 𝑠 1 ( 𝑡 ) = 1 − 𝑒 − ( 𝜔 0 𝑡 / 5 ) 2 s i n 2 𝜋 𝑓 0 𝑡 𝑠 , ( 5 a ) 2 ( 𝑡 ) = 1 − 𝑒 − ( 𝜔 0 𝑡 / 5 ) 2 s i n 2 𝜋 𝑓 0 𝑡 − 𝜋 / 2 , ( 5 b ) 𝑓 0 where 1 − 𝑒 − ( 𝜔 0 𝑡 / 5 ) 2 is the Larmor frequency (127.72 MHz for 1H at 3 T). The introduction of the coefficient ( 𝐵 1 ) in the excitation signal accelerated stabilization of the calculated 𝐸 𝑥 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 component 𝐸 𝑥 𝑦 was divided into 𝐸 𝑥 𝑧 and 𝐸 𝑥 = 𝐸 𝑥 𝑦 + 𝐸 𝑥 𝑧 , ( 6 ) : 𝐸 𝑛 + 1 𝑥 𝑦 1 𝑖 + 2 , 𝑗 , 𝑘 = 𝑒 − ( 𝜎 𝑦 / 𝜀 0 ) Δ 𝑡 𝐸 𝑛 𝑥 𝑦 1 𝑖 + 2 + 1 , 𝑗 , 𝑘 𝜎 𝑦 Δ 𝑦 1 − 𝑒 − ( 𝜎 𝑦 / 𝜀 0 ) Δ 𝑡 × 𝐻 𝑛 + 1 / 2 𝑧 𝑥 1 𝑖 + 2 1 , 𝑗 + 2 , 𝑘 − 𝐻 𝑛 + 1 / 2 𝑧 𝑥 1 𝑖 + 2 1 , 𝑗 − 2 − 𝐻 , 𝑘 𝑛 + 1 / 2 𝑧 𝑦 1 𝑖 + 2 1 , 𝑗 , 𝑘 + 2 − 𝐻 𝑛 + 1 / 2 𝑧 𝑦 1 𝑖 + 2 1 , 𝑗 , 𝑘 − 2 , 𝐸 ( 7 a ) 𝑛 + 1 𝑥 𝑧 1 𝑖 + 2 , 𝑗 , 𝑘 = 𝑒 − ( 𝜎 𝑧 / 𝜀 0 ) Δ 𝑡 𝐸 𝑛 𝑥 𝑧 1 𝑖 + 2 + 1 , 𝑗 , 𝑘 𝜎 𝑧 Δ 𝑧 1 − 𝑒 − ( 𝜎 𝑧 / 𝜀 0 ) Δ 𝑡 × 𝐻 𝑛 + 1 / 2 𝑦 𝑥 1 𝑖 + 2 1 , 𝑗 + 2 , 𝑘 − 𝐻 𝑛 + 1 / 2 𝑦 𝑥 1 𝑖 + 2 1 , 𝑗 − 2 − 𝐻 , 𝑘 𝑛 + 1 / 2 𝑦 𝑧 1 𝑖 + 2 1 , 𝑗 , 𝑘 + 2 − 𝐻 𝑛 + 1 / 2 𝑦 𝑧 1 𝑖 + 2 1 , 𝑗 , 𝑘 − 2 . ( 7 b ) where 𝐻 𝑥 Similarly, the magnetic field component 𝐻 𝑥 𝑦 was divided into 𝐻 𝑥 𝑧 and 𝐻 𝑥 = 𝐻 𝑥 𝑦 + 𝐻 𝑥 𝑧 , ( 8 ) : 𝐻 𝑛 + 1 / 2 𝑥 𝑦 1 𝑖 , 𝑗 + 2 1 , 𝑘 + 2 = 𝑒 − ( 𝜎 ∗ 𝑦 / 𝜇 0 ) Δ 𝑡 𝐻 𝑛 − 1 / 2 𝑥 𝑦 1 𝑖 , 𝑗 + 2 1 , 𝑘 + 2 + 1 𝜎 ∗ 𝑦 Δ 𝑦 1 − 𝑒 − ( 𝜎 ∗ 𝑦 / 𝜇 0 ) Δ 𝑡 × 𝐸 𝑛 𝑧 𝑥 1 𝑖 , 𝑗 + 1 , 𝑘 + 2 − 𝐸 𝑛 𝑧 𝑥 1 𝑖 , 𝑗 , 𝑘 + 2 − 𝐸 𝑛 𝑧 𝑦 1 𝑖 , 𝑗 + 1 , 𝑘 + 2 − 𝐸 𝑛 𝑧 𝑦 1 𝑖 , 𝑗 , 𝑘 + 2 , 𝐻 ( 9 a ) 𝑛 + 1 / 2 𝑥 𝑧 1 𝑖 , 𝑗 + 2 1 , 𝑘 + 2 = 𝑒 − ( 𝜎 ∗ 𝑦 / 𝜇 0 ) Δ 𝑡 𝐻 𝑛 − 1 / 2 𝑥 𝑧 1 𝑖 , 𝑗 + 2 1 , 𝑘 + 2 + 1 𝜎 ∗ 𝑦 Δ 𝑧 1 − 𝑒 − ( 𝜎 ∗ 𝑦 / 𝜇 0 ) Δ 𝑡 × 𝐸 𝑛 𝑦 𝑥 1 𝑖 , 𝑗 + 1 , 𝑘 + 2 − 𝐸 𝑛 𝑦 𝑥 1 𝑖 , 𝑗 , 𝑘 + 2 − 𝐸 𝑛 𝑦 𝑧 1 𝑖 , 𝑗 + 1 , 𝑘 + 2 ) − 𝐸 𝑛 𝑦 𝑧 1 𝑖 , 𝑗 , 𝑘 + 2 . ( 9 b ) where 𝑁 P M L 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, 𝑅 ( 0 ) , must be determined before the computations. Usually, given a required reflection coefficient at the outer layer with normal incidence (0° incidence angle), 𝜎 𝑖 = 𝑝 𝑖 𝑝 𝑚 𝜎 m a x , 𝜎 ∗ 𝑖 = 𝑝 𝑖 𝑝 𝑚 𝜎 ∗ m a x , ( 1 0 ) , σi, and σi* are obtained by [2, 3] 𝑝 where 𝑝 𝑖 is the thickness of the PML, 𝑖 is the distance from the inner side (towards the coil) of the boundary to 𝑚 th layer of PML, 𝜎 m a x is a coefficient in the range of 1 to 5, 𝜎 ∗ m a x and 𝜎 m a x = ( 𝑚 + 1 ) 𝜀 0 𝑐 0 1 2 𝑝 l n , 𝜎 𝑅 ( 0 ) ∗ m a x = 𝜀 0 𝜇 0 𝜎 m a x , ( 1 1 ) are the electric conductivity and magnetic conductivity of the outer layer of PML respectively, which are given by 𝜀 0 where 𝜇 0 and 𝑚 = 2 are the permittivity and permeability in free space, respectively.

According to Berenger [2, 3], an 𝑅 ( 0 ) = 1 0 − 4 and an 𝑝 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, “ 𝑠 ,” and the distance from the coil structure to the PML, “ 𝑝 ” (see Figure 2). We therefore use “ 𝑠 ” and “ 𝑁 𝑝 ” to assess the reflection of the PML. For quantitative measurements, we used “ 𝑝 ,” the layers of the PML, to denote the thickness of the PML, “ 𝑠 .” Each layer was modeled by one grid. The distance “ 𝑁 𝑠 ” was also measured in grids “ 𝑥 𝑜 𝑧 .”

fig2
Figure 2: The schematic view of the computational domain. (a) Longitudinal plane ( 𝑥 𝑜 𝑦 -plane in Figure 1(b)). (b) Transverse plane ( 𝐵 1 -plane in Figure 1(b)). “D” homogenous region of 𝑑 -field (shaded); “p”: the thickness of the PML; “s”: the distance from PML to the closest point on the coil geometry; “ 𝑙 ”: the diameter of the coil; “ 𝐵 1 ”: 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.

Usually, the 𝐷 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 “ 𝑑 ℎ ,” together with the diameter “ 𝑙 ℎ ” and the length “ 𝐵 1 ,” 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 calculated 𝑁 𝑐 = 𝑁 𝑚 , ( C 1 . a ) 𝑑 𝑓 𝑖 m a x = ( 𝑓 𝑐 𝑖 − 𝑓 𝑚 𝑖 ) m a x < 0 . 5 M H z , 𝑖 = 0 , 1 , … , 𝑁 𝑚 , ( C 1 . b ) 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 𝑁 𝑐 𝑁 𝑚 where 𝑑 𝑓 𝑖 is the number of the calculated modes; 𝑖 is the number of the measured modes; 𝑓 𝑐 𝑖 is the difference between the calculated and measured frequency at 𝑖 th mode; 𝑓 𝑚 𝑖 is the calculated resonance frequency at 𝑖 th mode; 𝐵 1 is the measured resonance frequency at 𝐵 1 th mode.

2. 𝐷 𝑟 m a x = 𝐵 1 𝑐 − 𝐵 1 𝑐 m a x 𝐵 1 𝑐 𝐵 < 1 0 % , 1 𝑐 ∈ 𝐷 , ( C 2 ) Field Homogeneity
We identified a region as homogenous region if the maximum relative deviation of the 𝐷 𝑟 m a x field there was <10%. This criterion, when applied to the birdcage coil, yielded 𝐵 1 𝑐 where 𝐵 1 is the maximum relative deviation; 𝐵 1 𝑐 is the calculated 𝐵 1 field; 𝐷 is the calculated mean 𝐵 1 field; 𝐵 1 is the cylindrical region (see Figure 2).

3. ( 𝐵 1 𝑐 𝑝 𝑗 − 𝐵 1 𝑐 𝑝 𝑘 ) m a x 𝐵 1 𝑐 𝑝 𝐵 < 1 0 % , 1 𝑐 ∈ 𝐷 ; 𝑁 𝑏 < 𝑗 , 𝑘 < 𝑁 𝑡 , ( C 3 ) Field Stability
The calculated 𝐵 1 𝑐 𝑝 𝑗 field should be stable after implementation of a certain number of time steps. We formulated the following criterion for this requirement: 𝐵 1 𝑐 𝑝 𝑘 where 𝑗 and 𝑘 are the 𝐵 1 th and 𝐵 1 𝑐 𝑝 th peak value of the calculated 𝐵 1 -field, respectively; 𝑁 𝑏 is the mean of the peak value of the calculated 𝐵 1 -field; 𝑁 𝑡 is the time steps at which the calculated 𝐵 1 -field starts to be stable; 𝑆 1 1 is the total time steps of calculation.

3. Results and Discussion

We found that both the resonance pattern and homogeneity of the calculated H 1 -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 Figure 3(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.

fig3
Figure 3: (a) Measured 𝑁 𝑝 = 5 of the coil (obtained with a network analyzer). 𝑁 𝑠 = 5 is the mode of interest (mode 1); (b) calculated resonance pattern with 𝑁 𝑝 = 1 0 and 𝑁 𝑠 = 8 grids. Only one mode resonated at 174 MHz; (c) calculated resonance pattern with 𝑆 1 1 and 𝑁 𝑝 = 1 0 grids. The number of modes and resonance frequencies is in excellent agreement with the measured 𝑁 𝑠 = 1 ; (d) calculated resonance pattern with H 1 and H 1 grids. Mode 0 and mode 1 are split (e.g., mode 1 split into two peaks: 𝑁 𝑝 and 𝑁 𝑠 .)

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, 𝑁 𝑝 and the distance from the PML to the closest surface of the coil, 𝑁 𝑠 , 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, 𝑁 𝑝 < 8 , must be greater than 8 and the number of grids between the closest point on the coil’s geometry and the PML interface, 𝑁 𝑠 < 5 , must be greater than 5. Any PML having an 𝐵 1 or an 𝐵 1 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. 𝐵 1 Homogeneity and Stability

To evaluate 𝑁 𝑝 homogeneity and stability, we first computed the 𝑁 𝑠 field of the coil using a PML with higher 𝑁 𝑝 and 𝑁 𝑠 values and then we repeated the computations using a PML with lower 𝑁 𝑝 and 𝑁 𝑠 values. Comparing the two computed results demonstrated that when the PML with higher 𝐵 1 and 𝐵 1 values increased the amplitude of the calculated 𝐵 1 -field gradually during the first 40 000 time steps (corresponding to about 40 cycles of the 𝐵 1 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 𝑁 𝑝 = 1 0 -field monotonically and without bound (see Figure 4(b)).

fig4
Figure 4: Calculated 𝑁 𝑠 = 1 0 -field at the center of the coil. (a) PML with 𝐵 1 and 𝑁 𝑝 = 5 . The 𝑁 𝑠 = 5 -field is stable after 40 000 time steps; (b) PML with 𝐵 1 and 𝐵 1 . The 𝑁 𝑝 -field increases monotonically even after 100 000 time steps.

These findings indicate that we obtained an optimal PML by iteratively calculating and evaluating the 𝑁 𝑠 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 𝑁 𝑝 equal 5, 8, 10, or 12 while increasing 𝑁 𝑝 from 1 to 20 for each value of 𝑁 𝑠 . While performing these numerical studies of the effects of varying 𝐵 1 and 𝐵 1 , we calculated the coil’s 𝑁 𝑝 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 𝑁 𝑠 field improved with increases of both 𝑁 𝑝 and 𝑁 𝑠 (see Figure 5). We noted that the calculated deviations were especially sensitive to values of 𝑁 𝑝 and 𝑁 𝑠 that were below 8 and 5, respectively, for 40-cycle computations, whereas increasing 𝐵 1 and 𝑁 𝑝 beyond 8 and 6, respectively, did not noticeably improve homogeneity of the calculated 𝑁 𝑠 field (see Figure 5(a)). When using 120 cycles, the minimum optimal 𝐵 1 and 𝑁 𝑝 increased to 10 layers and 12 grids (see Figure 5(b)).

fig5
Figure 5: Maximum deviations of the calculated 𝑁 𝑠 -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.

Because computational time is proportional to the number of grids, which in turn is directly related to the sum of 𝑁 𝑝 and 𝑁 𝑠 to the third power, 𝐵 1 and 𝐵 1 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 𝐵 1 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 𝐵 1 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 calculated 𝐵 1 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 𝐵 1 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.

fig6
Figure 6: Measurement of homogeneity of 𝑁 𝑝 = 1 0 field. (a) Acquired MR image of an oil phantom loaded within the coil. (b) Calculated 𝑁 𝑠 = 1 0 -map after 40-cycle iterations with 𝐵 1 and 𝐵 1 . The homogeneity of 𝑁 𝑝 = 1 0 -map is good inside, on the surface of, and outside of the coil. (c) Calculated 𝑁 𝑠 = 1 -map after 40-cycle iterations with 𝐵 1 and 𝐵 1 . The homogeneity of 𝑁 𝑝 = 1 0 -map is unacceptable. (d) Calculated 𝑁 𝑠 = 1 0 -map after 120-cycle iterations with 𝐵 1 and 𝐵 1 . The intensities of the 𝑄 field at the left-bottom corner and right-top corner are higher than those at other corners.
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 Table 1) [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 𝐵 1 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.

tab1
Table 1: The dielectric and conductivity of tissues inside a human head at 128 MHz.

The 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 calculated 𝑄 -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 𝐵 1 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 𝐵 1 -factor and deviation of the 𝐵 1 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 and 𝐵 1 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 𝐵 1 field inside the coil, and (3) the numerical stability of the calculated 𝐵 1 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 𝐵 1 field, a poor choice of PML parameters can produce a large fluctuation of the 𝐵 1 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 B 1 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 B 1 fields.

References

  1. K. S. Yee, “Numerical solution of initial boundary value problem involving Maxwell's equations in isotropic media,” IEEE Transaction on Antennas and Propagation, vol. 14, no. 3, pp. 302–307, 1966. View at Publisher · View at Google Scholar · View at Zentralblatt MATH
  2. J.-P. Berenger, “A perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 114, no. 2, pp. 185–200, 1994. View at Publisher · View at Google Scholar · View at MathSciNet
  3. J.-P. Berenger, “Three-dimensional perfectly matched layer for the absorption of electromagnetic waves,” Journal of Computational Physics, vol. 127, no. 2, pp. 363–379, 1996. View at Publisher · View at Google Scholar · View at MathSciNet
  4. Z. S. Sacks, D. M. Kingsland, R. Lee, and J.-F. Lee, “A perfectly matched anisotropic absorber for use as an absorbing boundary condition,” IEEE Transactions on Antennas and Propagation, vol. 43, no. 12, pp. 1460–1463, 1995. View at Publisher · View at Google Scholar
  5. S. D. Gedney, “An anisotropie perfectly matched layer-absorbing medium for the truncation of FDTD lattices,” IEEE Transactions on Antennas and Propagation, vol. 44, no. 12, pp. 1630–1639, 1996. View at Publisher · View at Google Scholar
  6. M. Kuzuoglu and R. Mittra, “Frequency dependence of the constitutive parameters of causal perfectly matched anisotropic absorbers,” IEEE Microwave and Guided Wave Letters, vol. 6, no. 12, pp. 447–449, 1996. View at Publisher · View at Google Scholar
  7. D. Correia and J.-M. Jin, “On the development of a higher-order PML,” IEEE Transactions on Antennas and Propagation, vol. 53, no. 12, pp. 4157–4163, 2005. View at Publisher · View at Google Scholar
  8. B. Sjögreen and N. A. Peterson, “Perfectly matched layers for Maxwell's equations in second order formulation,” Journal of Computational Physics, vol. 209, no. 1, pp. 19–46, 2005. View at Publisher · View at Google Scholar · View at MathSciNet
  9. J.-P. Bérenger, “Numerical reflection from FDTD-PMLS: a comparison of the split PML with the unsplit and CFS PMLS,” IEEE Transactions on Antennas and Propagation, vol. 50, no. 3, pp. 258–265, 2002. View at Publisher · View at Google Scholar
  10. J. Chen, Z. Feng, and J.-M. Jin, “Numerical simulation of SAR and B1-field inhomogeneity of shielded RF coils loaded with the human head,” IEEE Transactions on Biomedical Engineering, vol. 45, no. 5, pp. 650–659, 1998. View at Publisher · View at Google Scholar · View at PubMed
  11. T. S. Ibrahim, R. Lee, B. A. Baertlein, and P.-M. L. Robitaille, “B1 field homogeneity and SAR calculations for the birdcage coil,” Physics in Medicine and Biology, vol. 46, no. 2, pp. 609–619, 2001. View at Publisher · View at Google Scholar
  12. T. S. Ibrahim, R. Lee, B. A. Baertlein, Y. Yu, and P.-M. L. Robitaille, “Computational analysis of the high pass birdcage resonator: finite difference time domain simulations for high-field MRI,” Magnetic Resonance Imaging, vol. 18, no. 7, pp. 835–843, 2000. View at Publisher · View at Google Scholar
  13. U. D. Nguyen, J. S. Brown, I. A. Chang, J. Krycia, and M. S. Mirotznik, “Numerical evaluation of heating of the human head due to magnetic resonance imaging,” IEEE Transactions on Biomedical Engineering, vol. 51, no. 8, pp. 1301–1309, 2004. View at Publisher · View at Google Scholar · View at PubMed
  14. C. E. Hayes, W. A. Edelsterin, J. F. Schenck, O.M. Mueller, and M. Eash, “An efficient, highly homogeneous radiofrequency coil for whole-body imaging at 1.5T,” Journal of Magnetic Resonance, vol. 63, pp. 622–628, 1985.
  15. M. A. Stuchly and S. S. Stuchly, “Dielectric properties of biological substances—tabulated,” Journal of Microwave Power, vol. 15, no. 1, pp. 19–26, 1980.
  16. Y.-S. Tsuei, A. C. Cangellaris, and J. L. Prince, “Rigorous electromagnetic modeling of chip-to-package (first-level) interconnections,” IEEE Transactions on Components, Hybrids, and Manufacturing Technology, vol. 16, no. 8, pp. 876–883, 1993. View at Publisher · View at Google Scholar