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 phantom MRI images. We then computed the resonance pattern and 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 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) [1013]. 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 () 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 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 [1013]. 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 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 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 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).

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 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]: 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]: where is the maximum propagation velocity of the electromagnetic waves within the simulated region. In this case, , and , 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 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: where .

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 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 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 Field

After calculating the coil’s resonance pattern, we then computed the 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 where is the Larmor frequency (127.72 MHz for 1H at 3 T). The introduction of the coefficient () 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 : where Similarly, the magnetic field component was divided into and : where 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, , must be determined before the computations. Usually, given a required reflection coefficient at the outer layer with normal incidence (0° incidence angle), , σ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, is a coefficient in the range of 1 to 5, and are the electric conductivity and magnetic conductivity of the outer layer of PML respectively, which are given by where and are the permittivity and permeability in free space, respectively.

According to Berenger [2, 3], an 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 “.”

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 “,” 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 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; is the measured resonance frequency at th mode.

2. Field Homogeneity
We identified a region as homogenous region if the maximum relative deviation of the field there was <10%. This criterion, when applied to the birdcage coil, yielded where is the maximum relative deviation; is the calculated field; is the calculated mean field; is the cylindrical region (see Figure 2).

3. Field Stability
The calculated field should be stable after implementation of a certain number of time steps. We formulated the following criterion for this requirement: where and are the th and th peak value of the calculated -field, respectively; is the mean of the peak value of the calculated -field; is the time steps at which the calculated -field starts to be stable; is the total time steps of calculation.

3. Results and Discussion

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

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, , must be greater than 8 and the number of grids between the closest point on the coil’s geometry and the PML interface, , must be greater than 5. Any PML having an or an 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. 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 and values increased the amplitude of the calculated -field gradually during the first 40 000 time steps (corresponding to about 40 cycles of the 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 -field monotonically and without bound (see Figure 4(b)).

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 and , 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 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 and increased to 10 layers and 12 grids (see Figure 5(b)).

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, and 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 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 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 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 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.

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

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