Microwave heating of waterrich solvents is a widely used processing technique in research and applications. Highquality outcome requires a uniform temperature environment; which, in turn, depends on the balancing of a variety of effects taking place during the heating. Here, we show that two inherent effects, namely, polarizationcharge shielding and electromagnetic resonances, play a critical role in shaping up the field pattern in the heated water sample. Polarizationcharge shielding produces an internal electric field sensitive to the sample size, shape, and orientation. Internal electromagnetic resonances result in a widely varying electric field, while also allowing much deeper field penetration than the attenuation length to allow largescale treatment. The key to temperature uniformity, thus, lies in an optimized thermal flow to balance the nonuniform energy deposition. These complicated processes are examined in simulation and interpreted physically. It is shown that a spherical sample is most favorable for obtaining a high temperature uniformity mainly because of its rotational symmetry. This conclusion is significant in that prevailing sample vessels are mostly nonspherical.
I. INTRODUCTION
Water is one of the most abundant substances, both in existence and in use. Its properties are of general interest. In particular, water behavior under microwave heating is of importance to scientific research and applications. For example, with the advantages of volumetric heating and much shorter treatment time, microwaveassisted synthesis has grown to an active field of research,^{1–3} for which water is environmentally the cleanest solvent.^{4,5} As another example, microwave heating has been explored as a method for wastewater treatment in recent years.^{6} The dielectric properties of water have been well documented.^{7–9} There have also been a number of numerical and experimental studies on microwave heating of liquid water ranging from molecular dynamics^{10} to performance enhancement in a specific microwave structure.^{11–15} These studies reveal the complexity of the water heating processes. In particular, heat convection plays a key role in the resulting temperature distribution. As discussed below, more issues still merit further investigations. The current study, based on a simplified model, is complementary to earlier studies by focusing on the underlying physics and general trends.
Although microwave heating is volumetric in nature, temperature nonuniformity has been a major hurdle to numerous applications. The exposure to a nonuniform field is a clear cause for this difficulty. There are, however, some lessunderstood causes, which are inherent in nature in that they persist even in a perfectly uniform field. Two major causes are addressed in the current study: polarizationcharge shielding and electromagnetic resonances.
Exposed to an electromagnetic wave with Efield amplitude E_{0,} a dielectric object reacts with slight microscopic displacements of its bond electrons to form macroscopic polarization charges, which partially shield the incident E_{0}. As a result, the dielectric object's internal Efield can be significantly smaller than E_{0} (Ref. 16, Chap. 4). Furthermore, for a nonspherical object, the shielding effect can be highly dependent on its orientation with respect to E_{0},^{17,18} which constitutes a significant cause of the heating nonuniformity.
Independently, electromagnetic resonances can readily take place in a dielectric object of dimensions comparable to or greater than half of a wavelength in the dielectric. As an electromagnetic wave enters such a dielectric object, multiple reflections off the inside walls can constructively superpose into a resonant standing wave at discrete frequencies. While this can significantly boost the field strength, it also leads to a complicated field pattern and, therefore, a much greater field nonuniformity.
In Sec. II, we examine in detail the effects of field shielding and electromagnetic resonances in water in the presence of an incident plane wave at 2.45 GHz, a predominantly used frequency for microwave heating. The field pattern (hence the heating rate) is shown to be characterized by three distinct regimes of the sample size, while an elongated sample is subject to the additional complication of orientation. Section III incorporates the effects of thermal flows and compares the relative merits of spherical and elongated shapes over a range of sizes. The optimum shape is found to be spherical and two radius ranges are identified for stable operation. Realistic issues concerning water sphere heating are discussed in Sec. IV, followed by the summary, conclusion, and significance in Sec. V.
II. EFFECTS OF POLARIZATIONCHARGE SHIELDING AND ELECTROMAGNETIC RESONANCES ON THE HEATING RATE
A. Model and relevant parameters
In Eq. (3), h(x) is the ratio of the actual heating rate normalized to the heating rate under E_{0}. In a linear medium where ε is independent of E_{0}, E_{d}(x) is proportional to E_{0}. Thus, h(x) is independent of E_{0} and it gives the spatial profile of the heating rate for any E_{0}. Given E_{0}, ω, and ε′′, h(x) determines the actual heating rate through Eq. (2).
We now approximate the ε of waterrich solvent by the ε of pure water. For microwave heating, key factors to consider are the temperature dependence of ε (= ε′+iε′′), the wavelength in the dielectric (λ_{d}), and the power attenuation length (δ). The temperature (T) dependence of ε′/ε_{0} and ε′′/ε_{0} of water, given by Ref. 9, is plotted in Fig. 1(a) for the household microwave frequency of 2.45 GHz.
λ_{d} and δ are given by (Ref. 16, p. 314)
Based on the data in Fig. 1(a) and Eq. (6), we obtain λ_{d} and δ in water as functions of T [Fig. 1(b)] with the following features relevant to the current study.
Figure 1(b) shows a significant reduction of wavelength in water (λ_{d}) from λ_{free} (=12.25 cm) due to the relatively large ε′ of water [Fig. 1(a)]. Importantly, λ_{d} varies from λ_{d} = 1.39 cm at T = 25 °C to λ_{d} = 1.57 cm at T = 80 °C. The 11% increase will cause a significant upward shift of all resonant frequencies as T rises from 25 to 80 °C.
It can be seen in Fig. 1(b) that δ increases sharply with a higher T due to a much decreased ε′′ [Fig. 1(a)]. For example, δ = 1.8 cm at T = 25 °C and δ = 5.8 cm at and T = 80 °C. The 3.2fold increase in δ implies a much increased attenuation length at a higher T, hence beneficial to the heating of largesize objects.
B. Water sphere—three regimes in heating rate in terms of sphere radius
As the incident wave penetrates into the sphere, the fields inside are radically different with and without resonances. Assuming that the water sphere is centered at the origin of coordinates, the Mie solution is in the form of an infinite series. Based on the documented formulas in Ref. 21, Fig. 2 displays $ h \xaf$ [the mean heating rate, Eq. (4)] and σ_{h}/ $ h \xaf$ (the percentage spread of the heating rate) as functions of R for T = 25 and T = 80 °C. Figure 2 is characterized by three distinct regimes:

The quasistatic regime (R ≪ λ_{d}). In this regime, ω at 2.45 GHz is well below the lowest resonant frequency; hence, polarizationcharge shielding plays the dominant role. As R → 0, $ h \xaf$→1.4 × 10^{−3} at T = 25 °C and $ h \xaf$ → 2.2 × 10^{−3} at T = 80 °C with a negligible σ_{h}/ $ h \xaf$, indicating a muchreduced, but nearly uniform Efield inside the sphere due to polarizationcharge shielding. This is in good agreement with Eq. (7).
As R rises to values close to or above λ_{d}, electromagnetic resonances start to set in, each with a finite resonant width. Hence, at a given radius, the modes are more and less overlapping in the sphere. While there is no sharp boundary, we refer to those modes with >3 dB separation (i.e., a trough falling more than one half of the neighboring peaks) as strong resonances and those with a shallower trough as overlapping resonances. The former is dominated by a single mode with a strong sensitivity to input frequency fluctuations, while the latter by multiple modes with a reduced frequency sensitivity. The associated stability issue in regard to the input frequency fluctuations will be discussed at the end of Sec. III.

The sharpresonance regime (λ_{d}/2 < R < 2λ_{d}). The sphere interior is now dominated by sharp electromagnetic resonances, with the mode spectrum given by Fig. 2(a). It can be seen that, as T increases from 25 to 80 °C, the resonances become sharper due to the reduction of ε′′ [Fig. 1(a)]. More significantly, this is accompanied by a shift of the resonant peaks to a higher R due to the reduction of ε′. The amount of the shift (∼11%) is predicted by Fig. 1(a) and Eq. (6a). At a fixed R, the shift is large enough to cause a mode switching as T rises, suggesting that this is a relatively unstable regime during the heating.

The multipleresonance regime (R≫λ_{d}). As R increases to much beyond λ_{d}, the mode spectrum becomes so crowded that a few densely populated weak resonances overlap to result in a smaller $ h \xaf$ and smoother variations of $ h \xaf$ and σ_{h}/ $ h \xaf$ with R. This implies a return to stable operation.
Overall, the combined polarizationcharge shielding and resonance effects result in variations of $ h \xaf$ and σ_{h}/ $ h \xaf$ by orders of magnitude over the range of R under consideration. However, as will be shown in Sec. IV, this can be mostly compensated by optimized thermal flows and multiple wave paths.
Figure 2 indicates that, at the same frequency (2.45 GHz), the $ h \xaf$ and σ_{h}/ $ h \xaf$ profiles are highly dependent on the relative values of R and λ_{d}. This is also reflected in the Efield patterns under the plane wave in Eq. (1). The Efield patterns, representative of the three regimes, are plotted on the xz plane (Fig. 3) at 25 °C for three values of R. In a linear medium, the Efield amplitude can be normalized to E_{0}. The normalized value is given by the color code in logarithmic scale.
Figure 3(a) is in the quasistatic regime (R ≪ λ_{d}) showing a strong polarizationcharge shielding effect. In the sharpresonance regime, the first resonance at R = 0.68 cm is displayed in Fig. 3(b). Its field pattern (a pure TE_{10} mode) is in good agreement with the analytical prediction of Ref. 20. Figure 3(c) plots the field patterns for R = 6 cm in the multipleresonance regime. Although the incident field decreases exponentially toward the center because R ≫ δ (=1.8 cm), the resonant mode the incident field has excited on its path to the center can still build up to a significant field near the center, as shown in Fig. 3(c). This enables largescale heating.
Previous theoretical and experimental studies^{23} indicate that the scattering cross section of a plane wave by a water sphere is relatively weak in the Rayleigh limit [R ≪ λ_{free}, Figs. 3(a) and 3(b)], then increases sharply as R approaches λ_{free}. This explains the scattered field pattern outside the sphere in Fig. 3(c).
C. Water cylinder—additional complication due to cylinder orientation
Figures 2 and 3 apply to a water sphere. In reality, the heated sample is usually nonspherical in shape. In chemical synthesis, for example, the reaction vessel is typically cylindrical. Here we consider a plane wave incident onto a water cylinder of radius a and length L. This is a configuration not amenable to an analytical solution, so the 3D highfrequency structure simulator (HFSS, see Appendix, part 1) is used for numerical solutions.
Assume that the water cylinder, with a fixed lengthtoradius ratio of L/a = 10, is under the incident plane wave in Eq. (1). The cylindrical case (and nonspherical case in general,) is complicated by the orientation angle (θ) of the cylinder axis with the xdirected wave polarization. Figure 4 plots $ h \xaf$ and σ_{h}/ $ h \xaf$ at T = 25 °C (dashed lines) and T = 80 °C (solid lines) as functions of a for a cylinder lying on the x axis [θ = 0, Fig. 4(a)] and another one on the z axis [θ = 90°. Figure 4(b)]. Similar to the case of a water sphere, Fig. 4 displays three distinct regimes: the quasistatic regime (a ≪ λ_{d}), sharpresonance regime (λ_{d}/4 < a < λ_{d}) of highQ modes, and multipleresonance regime (a ≫ λ_{d}) of densely populated lowQ modes. Also, as T increases from 25 to 80 °C, there is a ∼11% shift of the resonant peaks to a higher R due to the reduction of ε′.
However, as a qualitative difference from the water sphere case, the water cylinder lacks the rotational symmetry. As a result, $ h \xaf$ is highly dependent on θ in the quasistatic regime, where $ h \xaf$(θ = 0) is ∼40 times greater than $ h \xaf$(θ = 90°). As discussed below, this is due to the strong orientational dependence of the polarizationcharge shielding effect (see Fig. 5, a = 0.1 cm data).
Figure 5 plots the overall Efield patterns in and around the water cylinder on the xz plane representative of the three regimes. Under the xpolarized, incident Efield, the polarization charges are separated in opposite Edirections. When the cylinder axis is aligned along the Efield (θ = 0), the ±polarization charges reside on opposite ends [e.g., Fig. 5(a)], hence farthest apart and producing the least shielding. When the cylinder axis is aligned perpendicular to E_{0} (θ = 90°), the ±polarization charges reside on opposite sidewalls [e.g., Fig. 5(b)], hence much closer and producing stronger shielding.
The shielding effect is most pronounced for small samples (e.g., Fig. 5, a = 0.1 cm). This is because, along with the incident Efield [Eq. (1)], the induced polarization charge varies sinusoidally in z with a period of λ_{f}. If the polarization charge is distributed over an area (at the end or on the sidewall) with a zdimension comparable to λ_{f} [e.g., Fig. 5(a), a = 3 cm] or greater than λ_{f} [e.g., Fig. 5(b), a = 3 cm], its sign variation along z will, thus, offset the overall shielding effect.
III. OPTIMUM SAMPLE SHAPE AND SIZE FOR TEMPERATURE UNIFORMITY
T(x) (hence $ T \xaf$ and σ_{T}) can be evaluated with the coupled software of HFSS and Icepak (described in Appendix), a suitable tool for the present purpose. Icepak returns the equilibrium solutions upon the specification of the geometrical configuration, the incident wave Efield amplitude E_{0} (or power density p_{inc} = $  E 0  2$/754 in SI unit), and the documented water properties listed in Table I. The initial temperature of both water and ambient air is assumed to be 25 °C.
Property .  Value .  Reference . 

Thermal conductivity (25 °C)  6.07 × 10^{−1} W/m K  Ref. 8 
Mass density (25 °C)  9.97 × 10^{2} kg/m^{3}  Ref. 8 
Specific heat (25 °C)  4.18 × 10^{3} J/kg K  Ref. 8 
Thermal expansion coefficient (volume)  2.56 × 10^{−4}/K  Ref. 6 
Radiation emissivity  9.50 × 10^{−1}  Ref. 22 
Temperaturedependent ε/ε_{0}  Fig. 1  Ref. 9 
Thermal diffusivity (25 °C)  1.46 × 10^{−7} m^{2}/s  Ref. 8 
Molecular mass  1.80 × 10^{−2 }kg/mol  Ref. 8 
Dynamic viscosity (25 °C)  8.90 × 10^{−4 }kg/m s  Ref. 8 
Property .  Value .  Reference . 

Thermal conductivity (25 °C)  6.07 × 10^{−1} W/m K  Ref. 8 
Mass density (25 °C)  9.97 × 10^{2} kg/m^{3}  Ref. 8 
Specific heat (25 °C)  4.18 × 10^{3} J/kg K  Ref. 8 
Thermal expansion coefficient (volume)  2.56 × 10^{−4}/K  Ref. 6 
Radiation emissivity  9.50 × 10^{−1}  Ref. 22 
Temperaturedependent ε/ε_{0}  Fig. 1  Ref. 9 
Thermal diffusivity (25 °C)  1.46 × 10^{−7} m^{2}/s  Ref. 8 
Molecular mass  1.80 × 10^{−2 }kg/mol  Ref. 8 
Dynamic viscosity (25 °C)  8.90 × 10^{−4 }kg/m s  Ref. 8 
A. Water sphere—effects of thermal conduction, convection, and radiation loss
We first consider the case of a water sphere by assuming a desired steadystate temperature of $ T \xaf$ = 80 °C for a certain application. This can be accomplished by adjusting the incident power (i.e., E_{0}). Figure 6 plots the required E_{0} to reach $ T \xaf$ = 80 °C and the resulting σ_{T}/ $ T \xaf$ as functions of the water sphere radius R. As expected, the required E_{0} is the greatest in the quasistatic regime (R ≪ λ_{d}) because of polarization charge shielding. It drops smoothly and significantly into the sharpresonance regime (λ_{d}/2 < R < 2λ_{d}), in which E_{0} fluctuates by a factor of ∼2 with R. The fluctuation is caused by a large change of the mode structure due to a small change in R. As R increases further into the multipleresonance regime (R ≫ λ_{d}), E_{0} is characterized again by a smooth variation with R.
Significantly, despite up to 100% heating rate spread [σ_{h}/ $ h \xaf$ in Fig. 2(b)], the steadystate temperature spread relative to the mean temperature [σ_{T}/ $ T \xaf$ in Fig. 6(b)] at $ T \xaf$ = 80 °C is rather low (0.2–1.1%) in all three regimes [Fig. 6(b)]. This is mainly due to heat convection and partly due to thermal conduction within the water sphere.
B. Water cylinder—orientational effect on the temperature spread
Similar to Fig. 6, we assume a desired steadystate temperature of $ T \xaf$ = 80 °C for a water cylinder with L/a = 10. Figure 7 plots the required E_{0} to reach $ T \xaf$ = 80 °C and the resulting σ_{T}/ $ T \xaf$ as functions of the cylinder radius a. For clarity, the simulated data points for a cylinder lying on the x axis (θ = 0) and z axis (θ = 90°) are connected by dashed and solid straight lines, respectively. The following differences are notable in comparison with the water sphere case (Fig. 6): (1) The results for a water cylinder rather strongly depend on the cylinder orientation. (2) The temperature spread (σ_{T}/ $ T \xaf$) is in general, greater than that for the water sphere.
C. Optimum sample shape and size
For a comparison of the characters and merits of spherical and cylindrical samples, portions of the data in Fig. 6 (spherical) and Fig. 7 (cylindrical) are superposed and plotted as functions of the reaction volume (V) in the range of 0.1–10 mL (Fig. 8). The range of V is chosen because most examples of microwaveassisted syntheses were performed on a less than 1 g scale (typical V = 1–5 mL).^{1} Figure 8 suggests a critical role of polarizationcharge shielding and indicates a significant advantage for the spherical geometry in terms of the temperature spread for reasons discussed below.

The Efield in cylindrical samples has a rather strong orientational dependence (Fig. 5). Furthermore, their sizes also significantly affect the temperature equilibration. As discussed in connection with Fig. 5, the shielding effect is most pronounced for small samples (e.g., Fig. 5, a = 0.1 cm). The surface area enclosing a unit volume is larger for smaller samples, especially near the end surface of the cylinder, where more heat (per unit volume) leaks into the air, hence resulting in a cooler spot. This is illustrated in the temperature profiles in Fig. 9 for a small (a = 0.1 cm) and a large (a = 0.5 cm) cylinder oriented along E_{0} (θ = 0), both at $ T \xaf$ = 80 °C. The smaller cylinder [Fig. 9(a)] has a centertoend temperature ratio of 1.20 while the larger cylinder has a smaller ratio of 1.08. Similarly (but not illustrated in figure), when the same two cylinders are oriented perpendicular to E_{0} (θ = 90°), the smaller one has a centertoend temperature ratio of 1.11 while the larger one has a smaller ratio of 1.01. This explains the rise of the temperature spread toward the smallV region in Fig. 8.

In comparison, for the (small) spherical samples in Fig. 8, the internal Efield is uniform [Eq. (7)]; hence, their equilibrium states feature a lower (or much lower) temperature spread than the cylindrical samples. Furthermore, the orientational independence of a sphere also makes the heating process much less sensitive to the applicator structure, hence more reproducible (discussed further in Sec. IV 1).
Focusing on the water sphere, there is still the issue of finding the optimum radius based on stability considerations. The magnetron is the predominantly used source for heating applications. As reported in Refs. 24 and 25, the measured spectrum of a household microwave oven is typically composed of a cluster of peaks spanning from 2.4 to 2.5 GHz. Furthermore, within this 4% band, it constantly changes with time and varies with the location, size, and geometry of the load.^{25} On the other hand, the required E_{0} for steadystate operation at a fixed $ T \xaf$ is sensitive to R in the sharpresonance regime (Fig. 6). A change in the operating frequency has the same effect as a proportional change in R. From Fig. 6, we find that a change of R up to 4% can present difficulties for equilibrium temperature control in the sharpresonance regime. This leaves two favorable radius ranges for stable operation: R < 0.7 cm (in the quasistatic regime) and R > 2 cm (in the multipleresonance regime). The percentage temperature spread is relatively small in both ranges (σ_{T}/ $ T \xaf$<1% for R < 0.7 cm and σ_{T}/ $ T \xaf$<0.3% for R > 2 cm). The R < 0.7 cm size (with a volume <1.2 mL) is suitable for research purpose, while the R > 2 cm size (with a volume >33.5 mL) is suitable for production purpose.
IV. FURTHER DISCUSSION ON A SPHERICAL SAMPLE

Applicability of the singlepath, travelingwave model to realistic situations
A great variety of 2.45 GHz microwave applicators (single or overmoded) are commercially available.^{26} The majority, if not all, employs a resonator structure, in which each standing wave is formed of two oppositely directed traveling waves bouncing back and forth between walls. Since a spherical sample is independent of the direction and polarization of the incident wave, the results obtained under the singlepath, travelingwave model (Fig. 6) apply to traveling waves from all directions in a resonator chamber. However, results for cylindrical samples under the singlepath, travelingwave model do not apply directly to an overmoded microwave cavity.

Deep field penetration through sphere resonances
The attenuation length (δ) of water at 25 °C is 1.8 cm [Fig. 1(b)]. Thus, the wave incident into a large sphere (R ≫ δ) damps sharply within a shallow outer region. However, densely populated resonant modes are excited as a result [Fig. 3(c)]. The natural structure of these normal modes extends to the whole sphere, hence allowing significant heating in the core region. Further, the attenuation length increases by a factor of 3.2 from 25 to 80 °C [Fig. 1(b)]. It is worth noting that, for these reasons, the percentage temperature spread (σ_{T}/ $ T \xaf$) is as low as 0.2% for R = 10 cm (≫δ) and shows a decreasing trend at still larger R [Fig. 6(b)].

Heating efficiency
Figure 6(a) shows that the incident E_{0} required for the steadystate operation at a fixed $ T \xaf$ varies widely with the sphere radius R. However, this does not preclude highefficiency heating over a wide range of R, because each traveling wave may execute multiple paths through the sphere in actual operations. Consider, for example, a water sphere being heated in a microwave resonator. The input structure (probe, loop, or aperture) can be tuned to the condition of critical coupling,^{27} under which all the incident power is fed into the cavity to be shared between the water sphere and the resonator wall. Since the wall is a good conductor, the water sphere normally absorbs the vast majority of the input energy.

The results reported here scale with R/λ_{d}. They can be generalized to other microwave frequencies or other solvents similar in dielectric properties.
V. SUMMARY, CONCLUSION, AND SIGNIFICANCE
Two basic physics effects, polarizationcharge shielding and electromagnetic resonances, are incorporated into the theory of microwave heating of water at the predominantly used frequency of 2.45 GHz. The former effect greatly reduces the internal Efield in water, while the latter effect significantly boosts the field strength along with a substantial field nonuniformity. A combination of these two effects results in ordersofmagnitude variations of the heating rate with respect to the sample size and shape (Figs. 2 and 4). However, this is only a snapshot of the overall process, which also involves thermal convection and conduction to even up the temperature. These thermal flows can lead to an equilibrium state with a high temperature uniformity.
A uniform temperature distribution requires measures to improve the heating uniformity. Unlike a sphere, an elongated object has a heating rate highly dependent on its orientation relative to the wave polarization (Fig. 4), which constitutes a major source for nonuniform heating. This is the principal reason for our conclusion that a spherical sample is most favorable for temperature uniformity (cf. Figs. 6 and 7). On the other hand, the incident wave can be muchattenuated in the core region of a large sphere. However, this can be compensated by the deep penetration of excited electromagnetic resonances [Fig. 4(c)].
Sharp electromagnetic resonances in the water sphere may, however, present a stability problem for a microwave source with poor frequency control (typically a magnetron). A quantitative examination has identified two radius ranges for stable operation: R < 0.7 and R > 2 cm (Sec. IV) suitable for research and largescale treatment, respectively.
The current study sheds light on key details, such as the thermal flow and internal temperature distribution, that are difficult to measure experimentally. The conclusion on the desirability of spherical sample for temperature uniformity is significant in that nonspherical (particularly cylindrical) reaction vessels are typically employed for laboratory tests. It is expected that the change to a spherical vessel may significantly improve the quality and reproducibility of the acquired data through a higher temperature uniformity.
As a final note, we have argued in Sec. III C about the relative advantages of a spherical sample as compared with a L/a = 10 cylinder. We have not studied other nonspherical samples. Since the sphere is unique in its orientational independence, these advantages are expected to carry over, to a greater or lesser extent, to other nonspherical shapes.
ACKNOWLEDGMENTS
The authors are grateful to ANSYS, Inc. for generously providing the HFSS and Icepak package as well as the technical support from Cybernet Systems Taiwan Co., Ltd. on its use. We are also grateful to Dr. H. H. Teng for helpful assistance. This work was funded by the National Science and Technology Council, Taiwan, under Grant No. MOST1112112M002016.
AUTHOR DECLARATIONS
Conflict of Interest
The authors have no conflicts to disclose.
Author Contributions
Li Chung Liu: Conceptualization (equal); Data curation (lead); Formal analysis (equal); Funding acquisition (supporting); Investigation (lead); Methodology (lead); Project administration (equal); Resources (lead); Software (lead); Supervision (equal); Validation (equal); Visualization (lead); Writing – original draft (equal); Writing – review & editing (equal). Ju Ching Liang: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (supporting); Investigation (equal); Methodology (equal); Project administration (supporting); Resources (equal); Software (lead); Supervision (supporting); Validation (equal); Visualization (equal); Writing – original draft (supporting); Writing – review & editing (equal). Kwan Wen Chen: Conceptualization (equal); Data curation (equal); Formal analysis (equal); Funding acquisition (supporting); Investigation (equal); Methodology (equal); Project administration (supporting); Resources (equal); Software (lead); Supervision (supporting); Validation (equal); Visualization (equal); Writing – original draft (supporting); Writing – review & editing (supporting). Kwo Ray Chu: Conceptualization (lead); Data curation (equal); Formal analysis (lead); Funding acquisition (lead); Investigation (lead); Methodology (lead); Project administration (lead); Resources (equal); Software (equal); Supervision (lead); Validation (lead); Visualization (lead); Writing – original draft (lead); Writing – review & editing (lead).
DATA AVAILABILITY
The data that support the findings of this study are available from the corresponding author upon reasonable request.
APPENDIX: SIMULATION SOFTWARES
1. Ansys HFSS
Electromagnetic field simulations were performed using the Ansys HFSS (highfrequency structure simulator) software, which solves the full set of Maxwell equations by the finite element method (http://anlage.umd.edu/HFSSv10UserGuide.pdf). In our simulations, PMLs (perfectly matched layers) are assumed on a closed boundary, which encloses the entire simulated structure. PMLs are fictitious materials that fully absorb the electromagnetic fields acting upon them.
2. Ansys ICEPAK
Thermal simulations were performed using the Ansys Icepak software. According to its manual (https://www.ansys.com/products/electronics/ansysicepak), the Icepak uses the computational fluid dynamics (CFD) solver engine for thermal and fluidflow calculations. Heat transfers through conduction, convection, and radiation. Laminar flow (assumed in the current study) is calculated by solving the Navier–Stokes equations for transport of mass, momentum, species, and energy to yield a fully coupled conduction and convection heat transfer prediction.
For the natural convection in a fluid which varies in temperature from one place to another, the gravity (through buoyant forces) drives a flow of fluid and heat transfer. In Icepak, the Boussinesq approximation is used to model the effect of density variations due to temperature while the fluid is assumed to be incompressible. The density is given by ρ = ρ_{0} − αρ_{0}(T − T_{0}), where ρ is the local density, ρ_{0} is the reference density, α is the coefficient of thermal expansion, T is the local temperature, and T_{0} is the reference temperature.
Radiation loss is calculated with the discrete ordinates (DO) radiation model, which solves the radiative transfer equation (RTE) for a finite number of discrete solid angles. Moreover, EM mapping and twoway coupling are used to perform twoway loss information between HFSS and Icepak. EM mapping allows Icepak to calculate the heat loss from electromagnetic fields in HFSS; twoway coupling allows HFSS to resimulate EM field with the varying temperature and Icepak to resimulate heat flow. HFSS and Icepak will do a coupled simulation for the steady state of either solid or fluid materials. The required input parameters are thermal conductivity, mass density, specific heat, thermal expansion coefficient, radiation emissivity, temperaturedependent electric properties (e.g., the dielectric constant) and, for fluid materials, the additional input parameters of thermal diffusivity, molecular mass, dynamic viscosity. Free openings are assumed on the boundary, which encloses the entire simulated structure, implying that the structure is exposed to free space.