RF engineers, oncologists, hyperthermia physicists.To compute the best achievable RF hyperthermia treatment performance using sources outside the head (“ultimate hyperthermia”) in realistic body models with non-uniform electrical and thermal properties and non-uniform tissue perfusion. Electromagnetic basis set: Our first step in computing the optimal external RF hyperthermia treatment is to obtain a basis set for electromagnetic (EM) fields in a realistic non-uniform body model. We have shown previously how to compute such a numerical basis [1] and used this approach to compute the ultimate signal-to-noise ratio in MRI [2]. First, we place a large number of dipoles (>105) around the body model at a distance no less than D=3 cm (Fig. 1). Second, we compute the incident E and H fields (i.e., in free space without the body model) created by random excitations of the dipole cloud. This is done by (i) exciting all the dipoles simultaneously using random excitations and (ii) computing the resulting E and H fields using the free-space Maxwell Green function. Finally, we compute the scattered E and H fields in the body model using the fast volume integral equation EM solver “MARIE” (https://github.com/thanospol/MARIE) [3]. Although bases of the solutions of Maxwell equations are strictly speaking infinite, we truncate the infinite basis to a finite number of basis vectors which converges to the ultimate with accuracy better than 1%. We model the Virtual Family “Duke” body model (head only, the body was cut off below the neck) [4] with a resolution of 3 mm isotropic (63×80×85 pixels), to which we added two tumors (one shallow, i.e. close to the skull, and one deep in the brain close to a ventricle). We compute EM basis sets in Duke at 64, 128, 297 and 600 MHz. Treatment optimization: We optimize the weights assigned to each basis vector (equivalent to voltage levels in a multi-channel hyperthermia coil) by maximization of the average specific absorption rate (SAR) inside the tumor while enforcing that the average SAR outside the tumor is below 3 W/kg. This is a convex constrained optimization problem that we solve using Matlab’s function “fmincon” using analytical expressions of the first and second derivatives of the objective and constraint functions for robust convergence. Temperature simulation: We solve Pennes bio-heat equation using the non-uniform thermal conductivity, heat capacity and perfusion maps of “Duke” using a Crank-Nicholson finite-difference time domain scheme with a time step computed using the Von Neuman stability criterion (code available at www.martinos.org/~guerin) [5]. Thermal parameter maps were obtained from the Gabriel database [6]. Perfusion maps were obtained from the Virtual Family database, except for tumors which were each decomposed in their necrotic core (perfusion 10% that of grey matter) and well-perfused penumbra (perfusion 200% that of grey matter) as shown in Fig. 3.Fig. 2 shows that convergence of the hyperthermia treatment efficacy (“SAR amplification ratio”) is achieved using ~500 basis vectors for all frequencies and both tumor depths. Roughly 230 and 100 basis vectors are needed to achieve 80% of the ultimate for the shallow and deep tumors, respectively. This suggests that, in practice, multi-channel hyperthermia coils with as many channels are needed to achieve results similar to the ultimate. The treatment frequency seems to affect the specificity of the ultimate treatment for the deep tumor only. This is likely because complex E fields interferences are required to “focus” the energy deposited in the deep tumor, which is more easily achieved at high frequency (=small wavelength). We point out that it is likely that the optimal treatment efficacy drops above a certain frequency however because of the shallower skin depth with increasing frequency. Fig. 3 shows SAR maps of the ultimate treatment at 600 MHz. Although the treatment energy is correctly focused in the tumor, there is considerable “spillover” to neighboring tissues. Note that we constrained the SAR outside the tumor to be below 3 W/kg, however the peak local SAR can – and obviously is – much greater in healthy tissues. For the deep tumor, the horn of the ventricle was significantly heated because of its proximity to the tumor and the high electrical conductivity of CSF. This suggests that an explicit SAR constraint corresponding to the ventricle volume should be added to the optimization algorithm. Temperature simulations in Fig. 4 reflect the SAR results of Fig. 3 and show a temperature amplification factor of x18 and x1.7 in the shallow and deep tumors, respectively, compared to the average healthy tissues.