Internal Temperature of Skin when Surface Temperature Is Controlled with an Electromagnetic Beam

Hongyun Wang¹, Wesley A. Burgei², Hong Zhou^3*

¹Department of Applied Mathematics University of California, Santa Cruz, USA.
²U.S. Department of Defense Joint Intermediate Force Capabilities Office, Quantico, USA.
³Department of Applied Mathematics Naval Postgraduate School, Monterey, USA.

DOI: 10.4236/jamp.2021.93030                

Abstract

We study the thermal effect on skin exposed to an electromagnetic beam of time-dependent power. We consider two types of beam power time schedules. In the controlled temperature exposure, the skin surface temperature is increased quickly to a prescribed level using a high beam power; then the surface temperature is maintained at the prescribed level by adjusting the beam power adaptively. In the constant power exposure, the applied beam power is relatively low and stays unchanged over the time. We start both types of exposures at the same time and compare their internal temperatures of skin when they have the same surface temperature. In a non-dimensionalized formulation, we show that at the moment when both exposure types reach the same prescribed surface temperature level, the controlled temperature exposure has a higher internal temperature at all depths. This conclusion is mathematically rigorous and is independent of skin material properties.

Keywords

Electromagnetic Beam of Time-Varying Power, Control of Surface Temperature, Thermal Effect of Beam, Non-Dimensional Analysis

1. Introduction

In many medical applications, such as cancer hyperthermia, patients are exposed to radiofrequency (RF) radiation [1]. We consider the thermal effect of RF radiation on human skin. The electromagnetic energy deposited by RF radiation increases the skin temperature. The skin surface temperature can be measured in real time by an infrared camera and can be controlled by tuning the power of RF radiation. In this study, we examine the skin internal temperature, which is the real target of control in medical applications.

We consider the situation where a test subject’s skin is exposed to an electromagnetic beam [2] [3] [4]. We focus on the thermal effect on the skin tissue along the beam center line, which has the highest power density deposited/absorbed relative to other lines parallel to the beam center line. On the skin surface and at any given depth, the beam center line has the highest temperature. In this mathematical study, we allow the beam power to vary with time. We consider two types of beam power time schedules. In the controlled temperature exposure, the skin surface temperature (at the beam center) is increased quickly to a prescribed level using a high beam power; then the surface temperature is maintained at the prescribed level by adjusting the beam power adaptively, leading to a time-varying beam power. The second type of exposure is the constant power exposure, in which the applied beam power is relatively low and stays unchanged over the time. We study the surface and internal temperatures of skin caused by these two types of exposures. We start both types of exposures at the same time. Due to its relatively higher initial power level, the controlled temperature exposure increases the temperature faster in the initial phase. Upon reaching the prescribed surface temperature level, the power is adaptively lowered to maintain the surface temperature. The constant power exposure, on the other hand, increases the temperature relatively slower but steadily. Eventually exposure types reach the prescribed surface temperature level. The main objective of this study is to compare the internal temperatures of the two exposures.

The rest of the article is organized as follows. In Section 2, we discuss the mathematical formulation and solution for the case where the beam power varies with time. Based on the solution for time-varying beam power, in Section 3 we develop the mathematical scheme for adjusting the beam power to maintain the surface temperature at the prescribed level. We run simulations to implement the control scheme and to demonstrate numerically that the beam power is a decreasing function of time in the controlled temperature exposure. This observation motivates Theorem 1. The internal temperatures of the two exposures are examined numerically in Section 4. A key observation is that when both exposures have the same surface temperature, the controlled temperature exposure always has the higher internal temperature at all depths. This finding motivates Theorem 2. In Section 5, we prove Theorem 1 and Theorem 2 rigorously in a dimensionless formulation. Thus, the main conclusions in this study are independent of skin material properties and independent of prescribed temperature level.

2. Mathematical Formulation for an Electromagnetic Beam of Time-Dependent Power

We adopt a formulation similar to the one in our previous studies [5] [6] [7], briefly summarized below. Let $T\left(z\mathrm{,}t\right)$ denote the skin temperature along the beam center line as a function of depth z and time t. We assume 1) the electromagnetic beam is perpendicular to the skin surface (i.e., beam incident angle = 0); 2) before the exposure, the skin has a uniform initial temperature $T_{\text{base}}$ (the baseline temperature); and 3) heat conduction is included only in the depth direction [8] (which is justified given the small length scale of electromagnetic wave penetrating in the skin depth direction and the much larger length scale of beam cross-section, and which allows us to separate z from $\left(x\mathrm{,}y\right)$ ).

The temperature distribution $T\left(z\mathrm{,}t\right)$ is governed by the heat equation

$\{\begin{array}{l}{\rho }_{\text{m}}{C}_p\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial z^2}+P\left(t\right)\mu \exp \left(-\mu z\right)\hfill \\ {\frac{\partial T}{\partial z}|}_{z=0}=\mathrm{0,}T\left(z\mathrm{,0}\right)=T_{\text{base}}\hfill \end{array}$ (1)

where

· ${\rho }_{\text{m}}$ is the mass density of the skin;

· $C_p$ is the specific heat capacity of the skin;

· $k$ is the heat conductivity of the skin;

· $\mu$ is the absorption coefficient of the skin for the beam frequency;

· $P\left(t\right)$ is the beam center power density absorbed into the skin at time t.

We first non-dimensionalize variables and functions in (1). The depth scale is provided by, $1/\mu$, which describes the characteristic scale of electromagnetic energy penetrating in the depth direction. The time scale is derived from the length scale and parameters of heat capacity and heat conduction. The temperature scale is usually set based on the objective of tests. For example, in studying heat-induced withdrawal reflex, the temperature scale is set to the difference between the activation temperature of nociceptors ( $T_{\text{act}}$ ) and the baseline temperature of skin ( $T_{\text{base}}$ ). The power density scale is derived from the temperature scale and skin material properties. Below we list these scales and the associated non-dimensional quantities.

· Depth scale and time scale:

$z_{\text{s}}\equiv \frac{1}{\mu }\mathrm{,}{t}_{\text{s}}\equiv \frac{{\rho }_{\text{m}}{C}_p}{k{\mu }^2}$

· Non-dimensional depth and time:

$z_{\text{nd}}\equiv \frac{z}{z_{\text{s}}}=\mu z\mathrm{,}{t}_{\text{nd}}\equiv \frac{t}{t_{\text{s}}}=t\frac{k{\mu }^2}{{\rho }_{\text{m}}{C}_p}$

· Temperature scale and power density scale:

$T_s\equiv T_{\text{act}}-T_{\text{base}}\mathrm{,}{P}_{\text{s}}\equiv k\mu T_s$

· Non-dimensional temperature as a function of $\left(z_{\text{nd}}\mathrm{,}{t}_{\text{nd}}\right)$ :

$T_{\text{nd}}\left(z_{\text{nd}}\mathrm{,}{t}_{\text{nd}}\right)\equiv \frac{T\left(z\mathrm{,}t\right)-T_{\text{base}}}{T_s}$

· Non-dimensional power density as a function of $t_{\text{nd}}$ :

$P_{\text{nd}}\left(t_{\text{nd}}\right)\equiv \frac{P\left(t\right)}{P_{\text{s}}}$

The non-dimensional temperature is governed by

$\{\begin{array}{l}\frac{\partial T_{\text{nd}}}{\partial t_{\text{nd}}}=\frac{{\partial }^2{T}_{\text{nd}}}{\partial z_{\text{nd}}^2}+P_{\text{nd}}\left(t_{\text{nd}}\right)\exp \left(-z_{\text{nd}}\right)\hfill \\ {\frac{\partial T_{\text{nd}}}{\partial z_{\text{nd}}}|}_{z_{\text{nd}}=0}=\mathrm{0,}T\left(z_{\text{nd}}\mathrm{,0}\right)=0\hfill \end{array}$ (2)

In this study, we analyze the non-dimensional system (2) and its solution. For conciseness, we shall drop the subscript “nd” and use the simple notations for all non-dimensional quantities. For example, $T\left(z\mathrm{,}t\right)$ means $T_{\text{nd}}\left(z_{\text{nd}}\mathrm{,}{t}_{\text{nd}}\right)$. The solution of initial boundary value problem (2) has the analytical expression

$\{\begin{array}{l}T\left(z\mathrm{,}t\right)={\displaystyle {\int }_0^t}P\left(s\right)G\left(z\mathrm{,}t-s\right)\text{d}s\hfill \\ G\left(z\mathrm{,}t\right)=\frac{1}{2}\text{erfc}\left(\frac{2t-z}{\sqrt{4t}}\right){\text{e}}^{t-z}+\frac{1}{2}\text{erfc}\left(\frac{2t+z}{\sqrt{4t}}\right){\text{e}}^{t+z}\hfill \end{array}$ (3)

where $\text{erfc}\left(u\right)$ is the complementary error function defined as

$\text{erfc}\left(u\right)\equiv \frac{2}{\sqrt{\pi }}{\displaystyle {\int }_u^{+\infty }}{\text{e}}^{-s^2}\text{d}s$ (4)

Next, we use the temperature solution (3) to design the beam power schedule $P\left(t\right)$ for controlling the skin surface temperature.

3. Surface Temperature Control

In this section, we study the controlled temperature exposure. Let $\Delta T$ be the prescribed surface temperature level. The $\Delta$ notation stems from that it is the intended non-dimensional surface temperature rise over the non-dimensional baseline temperature (0). Based on the temperature solution given in (3), we write the surface temperature as

$T\left(\mathrm{0,}t\right)={\displaystyle {\int }_0^t}P\left(s\right)b\left(t-s\right)\text{d}s$ (5)

where $b\left(t\right)\equiv G\left(\mathrm{0,}t\right)={\text{e}}^t\text{erfc}\left(\sqrt{t}\right)>0$

The asymptotic behavior of $b\left(t\right)$ follows from the asymptotic expansion of $\text{erfc}\left(\sqrt{t}\right)$.

$\text{erfc}\left(u\right)=\frac{{\text{e}}^{-u^2}}{\sqrt{\pi u^2}}\left(1-\frac{1}{2u^2}+\cdots \right)$

$b\left(t\right)={\text{e}}^t\text{erfc}\left(\sqrt{t}\right)=\frac{1}{\sqrt{\pi t}}\left(1-\frac{1}{2t}+\cdots \right)$ (6)

When the beam power is kept at any fixed value, $P\left(t\right)\equiv P_C$, the surface temperature is proportional to the beam power and increases monotonically with time without bound.

$T\left(\mathrm{0,}t\right)=P_C\beta (\; t\; )$

$\begin{array}{c}\beta \left(t\right)\equiv {\displaystyle {\int }_0^t}b\left(s\right)\text{d}s={\text{e}}^t\text{erfc}\left(\sqrt{t}\right)-1+\frac{2}{\sqrt{\pi }}\sqrt{t}\\ \approx \frac{2}{\sqrt{\pi }}\sqrt{t}-1+\frac{1}{\sqrt{\pi t}}\text{for}\text{large}t\mathrm{.}\end{array}$ (7)

In the controlled temperature exposure, we start with a relatively high beam power $P_0$. We keep the beam power at $P_0$ until time $t_0$ when the surface temperature reaches the prescribed $\Delta T$. Mathematically, $t_0$ is governed by $T\left(\mathrm{0,}{t}_0\right)=\Delta T$, which gives the equation

$\Delta T=P_0{\displaystyle {\int }_0^{t_0}}b\left(t_0-s\right)\text{d}s$ (8)

For $t>t_0$, the beam power is adjusted adaptively to maintain the surface temperature at the prescribed $\Delta T$. To mimic the realistic experimental situation, we consider the case of adjusting the beam power in discrete time steps. We use a uniform grid for $t>t_0$.

$t_j=t_0+j\Delta t,j=1,2,3,\cdots$

In each small time interval $\left(t_{j-1}\mathrm{,}{t}_j\right]$, we use a constant beam power $P_j$. The discrete-time beam power schedule has the form

$P\left(t\right)=\{\begin{array}{ll}{P}_0,\hfill & 0\le t\le t_0\hfill \\ P_j,\hfill & t\in \left(t_{j-1},t_j\right],j=1,2,3,\cdots \hfill \end{array}$ (9)

In the controlled temperature exposure, the initial beam power $P_0$ and the target temperature level $\Delta T$ are prescribed as the specified parameters. In comparison, the initial exposure period $t_0$ and the subsequent beam power levels $\left\{P_j,j=1,2,3,\cdots \right\}$ are the derived parameters: they are calculated from the given parameters. $t_0$ is determined from $T\left(0,t_0\right)=\Delta T$. For each $j\ge 1$, we set beam power $P_j$ to maintain $T\left(0,t_j\right)=\Delta T$. The beam power $P_j$ for time interval $\left(t_{j-1}\mathrm{,}{t}_j\right]$ depends on the power levels for all preceding time intervals. First, we solve for $P_1$ in equation $\Delta T=T\left(\mathrm{0,}{t}_1\right)$ :

$\Delta T=P_0{\displaystyle {\int }_0^{t_0}}b\left(t_1-s\right)\text{d}s+P_1{\displaystyle {\int }_{t_0}^{t_1}}b\left(t_1-s\right)\text{d}s$ (10)

Once $P_1$ is determined, we solve for $P_2$ in equation $\Delta T=T\left(\mathrm{0,}{t}_2\right)$ :

$\Delta T=P_0{\displaystyle {\int }_0^{t_0}}b\left(t_2-s\right)\text{d}s+P_1{\displaystyle {\int }_{t_0}^{t_1}}b\left(t_2-s\right)\text{d}s+P_2{\displaystyle {\int }_{t_1}^{t_2}}b\left(t_2-s\right)\text{d}s$ (11)

In general, when all preceding power levels $\left\{P_0\mathrm{,}{P}_1\mathrm{,}\cdots \mathrm{,}{P}_j\right\}$ are known, we solve for $P_{j+1}$ in equation $\Delta T=T\left(0,t_{j+1}\right)$ :

$\Delta T=P_0{\displaystyle {\int }_0^{t_0}}b\left(t_{j+1}-s\right)\text{d}s+{\displaystyle \underset{k=1}{\overset{j}{\sum }}}{P}_k{\displaystyle {\int }_{t_{k-1}}^{t_k}}b\left(t_{j+1}-s\right)\text{d}s+P_{j+1}{\displaystyle {\int }_{t_j}^{t_{j+1}}}b\left(t_{j+1}-s\right)\text{d}s$ (12)

It is worthwhile to compare the mathematical control described above and the feedback control implemented in real tests. In experiments, the skin material properties are unknown; the skin surface temperature is monitored in real time using an IR camera connected to a computer; and the computer switches the beam power on/off depending on the current surface temperature reading. At the beginning of each test, the high beam power $P_0$ is kept on until the prescribed temperature level $\Delta T$ is reached. Then beam power $P_0$ is switched on and off rapidly using the feedback control to maintain the surface temperature at $\Delta T$. When averaged over a time interval, the rapid on/off cycling of high beam power $P_0$ gives an effective constant (low) beam power over that time interval, which corresponds to beam power $P_j$ for time interval $\left(t_{j-1}\mathrm{,}{t}_j\right]$ in our mathematical formulation (9).

We first explore numerically the behavior of the beam power time schedule $P\left(t\right)$, constructed above for maintaining the surface temperature at $\Delta T$. Figure 1 plots $P\left(t\right)$ (left panel) and the associated surface temperature $T\left(\mathrm{0,}t\right)$ (right panel). The calculation is based on non-dimensional parameters $P_0=8$ and $\Delta T=2$. The initial exposure period is $t_0=0.3691$. All quantities are non-dimensional. The beam power drops to $P=3.1575$ at $t_0$ and from there it continues a gradual downward trend over the time. At $t=4.5$, the beam power is below $P=0.55$. Figure 1 suggests a key result regarding the beam power schedule calculated from our mathematical model for controlling surface temperature.

Theorem 1 Beam power levels $\left\{P_j\right\}$ calculated based on $\left\{T\left(0,t_j\right)=\Delta T\right\}$ satisfy

$P_0>P_1>P_2>\cdots >P_j>P_{j+1}>\cdots$ .

In other words,in the controlled temperature exposure,the beam power schedule $P\left(t\right)$ is a decreasing function of time.

This theorem is a key analytical tool when we compare the skin internal temperature for the two types of exposures: controlled temperature exposure vs constant beam power exposure. Although Theorem 1 is confirmed numerically in Figure 1 for $\Delta T=2$ and $P_0=8$, we will prove it rigorously in Section 5 for all values of $\Delta T$ and $P_0$. In Section 4, we compare the skin internal temperature of the two exposure types and summarize the key result in Theorem 2, which is also proved rigorously in Section 5.

Figure 1. The controlled temperature exposure. Left panel: beam power time schedule for maintaining surface temperature. Right panel: surface temperature vs time.

4. Skin Internal Temperature of Controlled Temperature Exposure vs Constant Power Exposure

In the constant power exposure [9], a relatively low beam power $P_L$ is applied over the time without any change in power level. To distinguish these two types of exposures, we use $P_{\Delta T}\left(t\right)$ and $P_L\left(t\right)$ to denote the beam power time schedule of respectively the controlled temperature exposure and the constant (low) power exposure.

$P_{\Delta T}\left(t\right)=\{\begin{array}{ll}{P}_0,\hfill & 0\le t\le t_0\hfill \\ P_j,\hfill & t\in \left(t_{j-1},t_j\right],j=1,2,3,\cdots \hfill \end{array}$ (13)

where $P_j$ is solved sequentially from

${\displaystyle {\int }_0^{t_{j-1}}}{P}_{\Delta T}\left(s\right)b\left(t_j-s\right)\text{d}s+P_j{\displaystyle {\int }_{t_{j-1}}^{t_j}}b\left(t_j-s\right)\text{d}s=\Delta T$

$P_L\left(t\right)\equiv P_L<P_0$, independent of t (14)

Accordingly, the corresponding skin internal temperature of these two exposure types are denoted respectively by $T_{\Delta T}\left(z\mathrm{,}t\right)$ and $T_L\left(z\mathrm{,}t\right)$.

$T_{\Delta T}\left(z\mathrm{,}t\right)={\displaystyle {\int }_0^t}{P}_{\Delta T}\left(s\right)G\left(z\mathrm{,}t-s\right)\text{d}s$

$T_L\left(z\mathrm{,}t\right)=P_L{\displaystyle {\int }_0^t}G\left(z\mathrm{,}t-s\right)\text{d}s$

Let $t_L$ be the time when the surface temperature of constant power exposure reaches the prescribed level. $t_L$ is governed by $T_L\left(\mathrm{0,}{t}_L\right)=\Delta T$, which via (7) becomes

$\Delta T=P_L\beta \left(t_L\right)$ (15)

$\beta \left(t_L\right)$ increases monotonically with $t_L$ without bound. It follows that for any $\Delta T>0$ and $P_L>0$, Equation (15) has a unique solution of $t_L$. When $\Delta T$ is fixed, $t_L$ is a decreasing function of $P_L$. Since $P_L<P_0$, we have $t_L>t_0$. That is, in the controlled temperature exposure, the surface temperature reaches $\Delta T$ earlier than it does in the constant power exposure. At time $t_L$, both the controlled temperature exposure and the constant power exposure have the same surface temperature ( $\Delta T$ ).

$T_{\Delta T}\left(0,t_L\right)=T_L\left(0,t_L\right)=\Delta T$ .

We are interested in comparing the skin internal temperature for the two exposure types at time $t_L$.

Again, we first explore it numerically. We use $P_0=8$ and $\Delta T=2$ for the controlled temperature exposure (the same parameters used in Figure 1). For the constant power exposure, we examine two power levels below $P_0$ : $P_L=0.8$ and $P_L=1.6$. The internal temperatures at time $t_L$ are shown in Figure 2 for the two exposure types and for the two values of $P_L$. Figure 2 suggests a key result regarding the internal temperatures of the two exposure types.

Figure 2. Comparison of internal temperatures $T_{\Delta T}\left(z\mathrm{,}{t}_L\right)$ (controlled temperature exposure) and $T_L\left(z\mathrm{,}{t}_L\right)$ (constant power exposure) at time $t_L$ when $T_L\left(\mathrm{0,}t\right)$ reaches $\Delta T$ and $T_{\Delta T}\left(\mathrm{0,}t\right)$ is maintained at $\Delta T$. Left panel: $P_L=0.8$. Right panel: $P_L=1.6$.

Theorem 2 When the surface temperature of the constant (low) power exposure reaches the prescribed,ΔT, both exposure types have the same surface temperature and the controlled temperature exposure always has a higher internal temperature.

$T_{\Delta T}\left(z,t_L\right)>T_L\left(z,t_L\right)\text{for}z>0$ .

We will prove Theorems 1 and 2 rigorously in the next section.

5. Proof of Theorems 1 and 2

An analytical expression of $T\left(z\mathrm{,}t\right)$ is given in (3). We rewrite it as

$\begin{array}{l}T\left(z\mathrm{,}t\right)={\displaystyle {\int }_0^t}P\left(s\right)b\left(t-s\right)g\left(z\mathrm{,}t-s\right)\text{d}s\\ b\left(t\right)\equiv G\left(\mathrm{0,}t\right)={\text{e}}^t\text{erfc}\left(\sqrt{t}\right)\mathrm{,}g\left(z\mathrm{,}t\right)\equiv \frac{G\left(z\mathrm{,}t\right)}{b\left(t\right)}\end{array}$ (16)

where $G\left(z\mathrm{,}t\right)$ is defined in (3). As mathematical preparation for the proof of Theorems 1 and 2, we study the properties of functions $G\left(z\mathrm{,}t\right)$, $b\left(t\right)$, and $g\left(z\mathrm{,}t\right)$.

5.1. Properties of $G\left(z,t\right)$, $b\left(t\right)$, and $g\left(z,t\right)$

Property 1 Function $G\left(z\mathrm{,}t\right)$ satisfies

$G\left(z,t\right)>G\left(0,t\right){\text{e}}^{\frac{-z^2}{4t}}\text{for}z>0$ .

Proof. We rewrite the integral in $\text{erfc}\left(u\right)$ using a change of variables $s=u+w$ :

$\text{erfc}\left(u\right)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_u^{+\infty }}{\text{e}}^{-s^2}\text{d}s=\frac{2}{\sqrt{\pi }}{\text{e}}^{-u^2}{\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-2uw-w^2}\text{d}w$ .

We apply this expression of $\text{erfc}(\cdot )$ to rewrite the two terms in (3).

$\begin{array}{c}\frac{1}{2}\text{erfc}\left(\frac{2t-z}{\sqrt{4t}}\right){\text{e}}^{t-z}=\frac{1}{\sqrt{\pi }}{\text{e}}^{t-z}{\text{e}}^{\frac{-{\left(2t-z\right)}^2}{4t}}{\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-\frac{2t-z}{\sqrt{t}}w-w^2}\text{d}w\\ =\frac{1}{\sqrt{\pi }}{\text{e}}^{\frac{-z^2}{4t}}{\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-\frac{2t-z}{\sqrt{t}}w-w^2}\text{d}w\end{array}$ .

$\frac{1}{2}\text{erfc}\left(\frac{2t+z}{\sqrt{4t}}\right){\text{e}}^{t+z}=\frac{1}{\sqrt{\pi }}{\text{e}}^{\frac{-z^2}{4t}}{\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-\frac{2t+z}{\sqrt{t}}w-w^2}\text{d}w$ .

Substituting these two terms back into (3), we get

$\begin{array}{c}G\left(z,t\right)=\frac{1}{\sqrt{\pi }}{\text{e}}^{\frac{-z^2}{4t}}{\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-\frac{2t}{\sqrt{t}}w-w^2}\underset{>2}{\underset{︸}{\left({\text{e}}^{\frac{z}{\sqrt{t}}w}+{\text{e}}^{-\frac{z}{\sqrt{t}}w}\right)}}\text{d}w\\ >\frac{1}{\sqrt{\pi }}{\text{e}}^{\frac{-z^2}{4t}}{\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-\frac{2t}{\sqrt{t}}w-w^2}2\text{d}w=G\left(0,t\right){\text{e}}^{\frac{-z^2}{4t}}\end{array}$ .

Here we have used $\left({\text{e}}^u+{\text{e}}^{-u}\right)>2$ for $u\ne 0$.

Property 2 $g\left(z\mathrm{,}t\right)$ is an increasing function of t at any fixed $z>0$.

Proof. We first calculate the time derivative of $G\left(z\mathrm{,}t\right)$. Differentiating (3) yields.

$\begin{array}{c}\frac{\partial G\left(z\mathrm{,}t\right)}{\partial t}=\frac{\partial }{\partial t}\left(\frac{1}{2}\text{erfc}\left(\frac{2t-z}{\sqrt{4t}}\right){\text{e}}^{t-z}+\frac{1}{2}\text{erfc}\left(\frac{2t+z}{\sqrt{4t}}\right){\text{e}}^{t+z}\right)\\ =G\left(z\mathrm{,}t\right)-\frac{1}{\sqrt{\pi }}{\text{e}}^{\frac{-z^2}{4t}}\frac{1}{\sqrt{t}}\end{array}$ (17)

Using (17) at z and at z = 0, we calculate the time derivative of $g\left(z\mathrm{,}t\right)=\frac{G\left(z\mathrm{,}t\right)}{G\left(\mathrm{0,}t\right)}$.

$\begin{array}{c}\frac{\partial g\left(z\mathrm{,}t\right)}{\partial t}=\frac{\partial }{\partial t}\left(\frac{G\left(z\mathrm{,}t\right)}{G\left(\mathrm{0,}t\right)}\right)=\frac{G\left(\mathrm{0,}t\right)\frac{\partial G\left(z\mathrm{,}t\right)}{\partial t}-G\left(z\mathrm{,}t\right)\frac{\partial G\left(\mathrm{0,}t\right)}{\partial t}}{G{\left(\mathrm{0,}t\right)}^2}\\ =\frac{\frac{1}{\sqrt{\pi }\sqrt{t}}\underset{>0}{\underset{︸}{\left(-G\left(0,t\right){\text{e}}^{\frac{-z^2}{4t}}+G\left(z\mathrm{,}t\right)\right)}}}{G{\left(\mathrm{0,}t\right)}^2}>0\end{array}$ .

Here we have used Property 1 to conclude that the term over the underbrace is positive.

Property 3 Let $f\left(t\right)\equiv \frac{b\left(t-a\right)}{b\left(t\right)},a>0$. Then $f\left(t\right)$ is a decreasing function of t for $t\ge a$.

Proof. We first show several items about function $b\left(t\right)$.

1) $b\left(t\right)={\text{e}}^t\text{erfc}\left(\sqrt{t}\right)=\frac{2}{\sqrt{\pi }}{\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-2\sqrt{t}w-w^2}\text{d}w$ decreases with t.

2) $\sqrt{t}b\left(t\right)=\frac{2\sqrt{t}}{\sqrt{\pi }}{\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-2\sqrt{t}w-w^2}\text{d}w=\frac{1}{\sqrt{\pi }}{\displaystyle {\int }_0^{+\infty }}{\text{e}}^{-u-\frac{u^2}{4t}}\text{d}u$ increases with t.

3) $b^{\prime }\left(t\right)=\left({\text{e}}^t\text{erfc}\left(\sqrt{t}\right)\right)=b\left(t\right)-\frac{1}{\sqrt{\pi }\sqrt{t}}$.

Graphs of $b\left(t\right)$ and $\sqrt{t}b\left(t\right)$ are illustrated in the left panel of Figure 3. The right panel of Figure 3 compares $g\left(z\mathrm{,}t\right)$ vs z for several values of t (Property 2). With the results itemized above, we examine the derivative of $f\left(t\right)$.

$\begin{array}{c}{f}^{\prime }\left(t\right)=\frac{b^{\prime }\left(t-a\right)b\left(t\right)-b\left(t-a\right)b^{\prime }\left(t\right)}{b{\left(t\right)}^2}\\ =\frac{\frac{1}{\sqrt{\pi }}\left(-\frac{1}{\sqrt{t-a}}b\left(t\right)+b\left(t-a\right)\frac{1}{\sqrt{t}}\right)}{b{\left(t\right)}^2}\\ =\frac{\sqrt{\left(t-a\right)}b\left(t-a\right)-\sqrt{t}b\left(t\right)}{\sqrt{\pi }\sqrt{t-a}\sqrt{t}b{\left(t\right)}^2}<0\end{array}$

Here we have used items 3 and 2 above.

Property 4 Suppose $h_1\left(t\right)$ and $h_2\left(t\right)$ are positive and decreasing functions of t, Let

$f\left(t\right)\equiv \frac{{\displaystyle {\int }_t^{t+\Delta t}}{h}_1\left(s\right)h_2\left(s\right)\text{d}s}{{\displaystyle {\int }_t^{t+\Delta t}}{h}_2\left(s\right)\text{d}s},\Delta t>0$

Then $f\left(t\right)$ is a decreasing function of t.

Proof. We examine the derivative of $f\left(t\right)$.

$f^{\prime }\left(t\right)=\frac{-\left(h_1\left(t\right)h_2\left(t\right)-h_1\left(t+\Delta t\right)h_2\left(t+\Delta t\right)\right){\displaystyle {\int }_t^{t+\Delta t}}{h}_2\left(s\right)\text{d}s+\left(h_2\left(t\right)-h_2\left(t+\Delta t\right)\right){\displaystyle {\int }_t^{t+\Delta t}}{h}_1\left(s\right)h_2\left(s\right)\text{d}s}{{\left({\displaystyle {\int }_t^{t+\Delta t}}{h}_2\left(s\right)\text{d}s\right)}^2}$.

Since both $h_1\left(t\right)$ and $h_2\left(t\right)$ are positive and decreasing, we have

${\displaystyle {\int }_t^{t+\Delta t}}{h}_1\left(s\right)h_2\left(s\right)\text{d}s<h_1\left(t\right){\displaystyle {\int }_t^{t+\Delta t}}{h}_2\left(s\right)\text{d}s$

$h_2\left(t\right)-h_2\left(t+\Delta t\right)>0$

$\begin{array}{l}{h}_1\left(t\right)h_2\left(t\right)-h_1\left(t+\Delta t\right)h_2\left(t+\Delta t\right)\\ =h_1\left(t\right)\left(h_2\left(t\right)-h_2\left(t+\Delta t\right)\right)+\left(h_1\left(t\right)-h_1\left(t+\Delta t\right)\right)h_2\left(t+\Delta t\right)\\ >h_1\left(t\right)\left(h_2\left(t\right)-h_2\left(t+\Delta t\right)\right)\end{array}$.

It follows that

$f^{\prime }\left(t\right)<\frac{\left(h_2\left(t\right)-h_2\left(t+\Delta t\right)\right)\left(-h_1\left(t\right){\displaystyle {\int }_t^{t+\Delta t}}{h}_2\left(s\right)\text{d}s+{\displaystyle {\int }_t^{t+\Delta t}}{h}_1\left(s\right)h_2\left(s\right)\text{d}s\right)}{{\left({\displaystyle {\int }_t^{t+\Delta t}}{h}_2\left(s\right)\text{d}s\right)}^2}<0$.

Property 5 Let

$f\left(t\right)\equiv \frac{{\displaystyle {\int }_0^{\Delta t}}b\left(t-a+s\right)\text{d}s}{{\displaystyle {\int }_0^{\Delta t}}b\left(t+s\right)\text{d}s},\Delta t>0,a>0$.

Then $f\left(t\right)$ is a decreasing function of t for $t\ge a$.

Figure 3. Left panel: graphs of $b\left(t\right)$ and $\sqrt{t}b\left(t\right)$. Right panel: $g\left(z\mathrm{,}t\right)$ vs z.

Proof. We cast $f\left(t\right)$ into the form of Property 4

$f\left(t\right)=\frac{{\displaystyle {\int }_t^{t+\Delta t}}{h}_1\left(s\right)h_2\left(s\right)\text{d}s}{{\displaystyle {\int }_t^{t+\Delta t}}{h}_2\left(s\right)\text{d}s}$

$h_2\left(s\right)=b\left(s\right),h_1\left(s\right)=\frac{b\left(s-a\right)}{b\left(s\right)}$.

Property 3 gives that both $h_1\left(s\right)$ and $h_2\left(s\right)$ are positive and decreasing functions of s. Consequently the decreasing of $f\left(t\right)$ follows directly from Property 4.

Property 6 Let

$\varphi \left(t\mathrm{<mo>\; ;\; </mo>}a\right)\equiv \frac{{\displaystyle {\int }_a^{a+\Delta t}}b\left(t-s\right)\text{d}s-{\displaystyle {\int }_0^{\Delta t}}b\left(t-s\right)\text{d}s}{{\displaystyle {\int }_a^{a+\Delta t}}b\left(t-s\right)\text{d}s}\mathrm{,}a>0,\Delta t>0$ (18)

Then $\varphi \left(t;a\right)<1$ and $\varphi \left(t\mathrm{<mo>\; ;\; </mo>}a\right)$ is a decreasing function of t for $t\ge \left(a+\Delta t\right)$.

Proof. We cast $\varphi \left(t\mathrm{<mo>\; ;\; </mo>}a\right)$ into the form of Property 5

$\varphi \left(t;a\right)=1-\frac{1}{f\left(t\right)},f\left(t\right)=\frac{{\displaystyle {\int }_0^{\Delta t}}b\left(t-\left(a+\Delta t\right)+s\right)\text{d}s}{{\displaystyle {\int }_0^{\Delta t}}b\left(t-\Delta t+s\right)\text{d}s}$

$\varphi \left(t;a\right)<1$ follows from that $f\left(t\right)$ is positive. Property 5 tells us that $f\left(t\right)$ decreases with t, which implies that $\varphi \left(t\mathrm{<mo>\; ;\; </mo>}a\right)$ is a decreasing function of t.

5.2. Proof of Theorem 1

We need to show $P_0<P_1<P_2<\cdots$. Since our main focus is on the relation among $\left\{P_j\right\}$, we like to write $P_{j+1}$ in terms of the preceding power levels. For $P_1$, we take the difference between (10) and (8) to obtain

$\begin{array}{c}{P}_1{\displaystyle {\int }_{t_0}^{t_1}}b\left(t_1-s\right)\text{d}s=P_0\left({\displaystyle {\int }_0^{t_0}}b\left(t_0-s\right)\text{d}s-{\displaystyle {\int }_0^{t_0}}b\left(t_1-s\right)\text{d}s\right)\\ =P_0\left({\displaystyle {\int }_{t_0}^{t_0+\Delta t}}b\left(t_1-s\right)\text{d}s-{\displaystyle {\int }_0^{\Delta t}}b\left(t_1-s\right)\text{d}s\right)\end{array}$.

It follows that

$P_1=P_0\varphi \left(t_1;t_0\right)<P_0$ (19)

where $\varphi \left(t\mathrm{<mo>\; ;\; </mo>}a\right)$ is defined in (18). For $P_2$, we take the difference between (11) and (10).

$\begin{array}{c}{P}_2{\displaystyle {\int }_{t_1}^{t_2}}b\left(t_2-s\right)\text{d}s=P_1\left({\displaystyle {\int }_{t_0}^{t_1}}b\left(t_1-s\right)\text{d}s-{\displaystyle {\int }_{t_0}^{t_1}}b\left(t_2-s\right)\text{d}s\right)\\ +P_0\left({\displaystyle {\int }_0^{t_0}}b\left(t_1-s\right)\text{d}s-{\displaystyle {\int }_0^{t_0}}b\left(t_2-s\right)\text{d}s\right)\end{array}$

which leads to

$\begin{array}{c}{P}_2=P_1\varphi \left(2\Delta t;\Delta t\right)+\left(1-\varphi \left(2\Delta t;\Delta t\right)\right)P_0\varphi \left(t_2;t_0\right)\\ =P_1+\left(1-\varphi \left(2\Delta t;\Delta t\right)\right)\underset{<0}{\underset{︸}{\left(-P_1+P_0\varphi \left(t_2;t_0\right)\right)}}<P_1\end{array}$ (20)

Here we have used Property 6 and Equation (19) to conclude $P_0\varphi \left(t_2;t_0\right)<P_0\varphi \left(t_1;t_0\right)=P_1$. For P₃, we take the difference between (12) and (11),

$\begin{array}{c}{P}_3{\displaystyle {\int }_{t_2}^{t_3}}b\left(t_3-s\right)\text{d}s=P_2\left({\displaystyle {\int }_{t_1}^{t_2}}b\left(t_2-s\right)\text{d}s-{\displaystyle {\int }_{t_1}^{t_2}}b\left(t_3-s\right)\text{d}s\right)\\ +P_1\left({\displaystyle {\int }_{t_0}^{t_1}}b\left(t_2-s\right)\text{d}s-{\displaystyle {\int }_{t_0}^{t_1}}b\left(t_3-s\right)\text{d}s\right)\\ +P_0\left({\displaystyle {\int }_0^{t_0}}b\left(t_2-s\right)\text{d}s-{\displaystyle {\int }_0^{t_0}}b\left(t_3-s\right)\text{d}s\right)\end{array}$

We introduce short notation ${\varphi }_k\equiv \varphi \left(k\Delta t\mathrm{<mo>\; ;\; </mo>}\Delta t\right)$ and write $P_3$ as

$\begin{array}{l}{P}_3=P_2{\varphi }_2+\left(1-{\varphi }_2\right)P_1{\varphi }_3+\left(1-{\varphi }_2\right)\left(1-{\varphi }_3\right))P_0\varphi \left(t_3;t_0\right)\\ =P_2+\left(1-{\varphi }_2\right)\underset{<0}{\underset{︸}{\left(-P_2+P_1{\varphi }_2+\left(1-{\varphi }_2\right)P_0\varphi \left(t_3;t_0\right)\right)}}\\ +\left(1-{\varphi }_2\right)\left({\varphi }_2-{\varphi }_3\right)\underset{<0}{\underset{︸}{\left(-P_1+P_0\varphi \left(t_3;t_0\right)\right)}}<P_2\end{array}$ (21)

Here we have used $\left({\varphi }_2-{\varphi }_3\right)>0$ and $\varphi \left(t_3;t_0\right)<\varphi \left(t_2;t_0\right)<\varphi \left(t_1;t_0\right)$ from Property 6. Continuing in this way, we can show $P_{j+1}<P_j$ by induction. In summary, in the controlled temperature exposure, the beam power time schedule is a decreasing function of time.

5.3. Proof of Theorem 2

We need to show $T_{\Delta T}\left(z,t_L\right)>T_L\left(z,t_L\right)$ for $z>0$. Here the time instance $t_L$ is defined by $T_L\left(z,t_L\right)=\Delta T$, which leads to

$0=T_{\Delta T}\left(0,t_L\right)-T_L\left(0,t_L\right)={\displaystyle {\int }_0^{t_L}}\left(P_{\Delta T}\left(s\right)-P_L\right)b\left(t_L-s\right)\text{d}s$ (22)

From Theorem 1, $P_{\Delta T}\left(t\right)$, the beam power schedule of controlled temperature exposure, is a decreasing function of t. Since $b\left(t_L-s\right)$ is positive, for the integral in (22) to be zero, there must be a unique $t^{\mathrm{<mo>\; *\; </mo>}}\in \left(\mathrm{0,}{t}_L\right)$ such that

$\{\begin{array}{l}{P}_{\Delta T}\left(t\right)>P_L\text{for}t<t^{*}\\ P_{\Delta T}\left(t\right)<P_L\text{for}t>t^{*}\end{array}$ (23)

Equation (22) is based on the difference in the surface temperature ( $z=0$ ) between the two exposure types. For the difference in the internal temperature of skin ( $z>0$ ), we use the expression of $T\left(z\mathrm{,}t\right)$ given in terms of $b\left(t\right)g\left(z\mathrm{,}t\right)$ in (16).

$\begin{array}{c}{T}_{\Delta T}\left(z,t_L\right)-T_L\left(z,t_L\right)={\displaystyle {\int }_0^{t^{*}}}\underset{>0}{\underset{︸}{\left(P_{\Delta T}\left(s\right)-P_L\right)}}b\left(t_L-s\right)\underset{>g\left(z,t_L-t^{*}\right)}{\underset{︸}{g\left(z,t_L-s\right)}}\text{d}s\\ +{\displaystyle {\int }_{t^{*}}^{t_L}}\underset{<0}{\underset{︸}{\left(P_{\Delta T}\left(s\right)-P_L\right)}}b\left(t_L-s\right)\underset{<g\left(z,t_L-t^{*}\right)}{\underset{︸}{g\left(z,t_L-s\right)}}\text{d}s\end{array}$ (24)

Property 2 tells us that $g\left(z\mathrm{,}t\right)$ is an increasing function of t.

$\{\begin{array}{l}g\left(z,t_L-s\right)>g\left(z,t_L-t^{*}\right)\text{for}s<t^{*}\\ g\left(z,t_L-s\right)<g\left(z,t_L-t^{*}\right)\text{for}s>t^{*}\end{array}$ (25)

Substituting (23) and (25) into (24), we obtain

$\begin{array}{l}{T}_{\Delta T}\left(z,t_L\right)-T_L\left(z,t_L\right)\\ >g\left(z,t_L-t^{*}\right){\displaystyle {\int }_0^{t^{*}}}\left(P_{\Delta T}\left(s\right)-P_L\right)b\left(t_L-s\right)\text{d}s\\ +g\left(z,t_L-t^{*}\right){\displaystyle {\int }_{t^{*}}^{t_L}}\left(P_{\Delta T}\left(s\right)-P_L\right)b\left(t_L-s\right)\text{d}s\\ =g\left(z,t_L-t^{*}\right)\underset{=T_{\Delta T}\left(0,t_L\right)-T_L\left(0,t_L\right)}{\underset{︸}{{\displaystyle {\int }_0^{t_L}}\left(P_{\Delta T}\left(s\right)-P_L\right)b\left(t_L-s\right)\text{d}s}}=0\end{array}$ (26)

In conclusion, when both exposure types reach the same surface temperature, the controlled temperature exposure always has a higher internal temperature.

6. Concluding Remarks

In this study, we considered the thermal effect on skin exposed to an electromagnetic beam. We investigated the skin surface temperature and internal temperature caused by the beam. Specifically, two exposure types were examined. In the controlled temperature exposure, a high beam power is used to increase the skin surface temperature quickly to a prescribed level. Then the beam power is adjusted adaptively to maintain the surface temperature at the prescribed level. In the constant power exposure, a relatively low beam power is applied without any change in power level over the time. We start both types of exposures at the same time. The controlled temperature exposure will reach the prescribed surface temperature level first since it has a higher initial beam power. To maintain the surface temperature once the prescribed level is attained, beam power drops significantly and keeps declining gradually over the time. When both types of exposures reach the same surface temperature, the controlled temperature exposure always has a higher internal temperature at all depths of skin. We proved this conclusion rigorously in a dimensionless formulation. This conclusion is independent of skin material properties, initial beam power levels and the prescribed surface temperature level.

Acknowledgement and Disclaimer

The authors acknowledge the Joint Intermediate Force Capabilities Office of U.S. Department of Defense and the Naval Postgraduate School for supporting this work. The views expressed in this document are those of the authors and do not reflect the official policy or position of the Department of Defense or the U.S. Government.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Ryan, K.L., D’Andrea, J.A. Jauchem, J.R. and Mason, P.A. (2000) Radio Frequency Radiation of Millimeter Wave Length: Potential Occupational Safety Issues Relating to Surface Heating. Health Physics, 78, 170-181. https://doi.org/10.1097/00004032-200002000-00006
[2]	Walters, T.J., Blick, D.W., Johnson, L.R., Adair, E.R. and Foster, K.R. (2000) Heating and Pain Sensation Produced in Human Skin by Millimeter Waves: Comparison to a Simple Thermal Model. Health Physics, 78, 259-267. https://doi.org/10.1097/00004032-200003000-00003
[3]	Zhadobov, M., Chahat, N., Sauleau, R., Le Quement, C. and Le Drean, Y. (2011) Millimeter-Wave Interactions with the Human Body: State of Knowledge and Recent Advances. International Journal of Microwave and Wireless Technologies, 3, 237-247. https://doi.org/10.1017/S1759078711000122
[4]	Topfer, F. and Oberhammer, J. (2015) Millimeter-Wave Tissue Diagnosis: The most Promising Fields for Medical Applications. IEEE Microwave Magazine, 16, 97-113. https://doi.org/10.1109/MMM.2015.2394020
[5]	Wang, H., Burgei, W.A. and Zhou, H. (2020) A Concise Model and Analysis for Heat-Induced Withdrawal Reflex Caused by Millimeter Wave Radiation. American Journal of Operations Research, 10, 31-81. https://doi.org/10.4236/ajor.2020.102004
[6]	Wang, H., Burgei, W.A. and Zhou, H. (2020) Non-Dimensional Analysis of Thermal Effect on Skin Exposure to an Electromagnetic Beam. American Journal of Operations Research, 10, 147-162. https://doi.org/10.4236/ajor.2020.105011
[7]	Wang, H., Burgei, W.A. and Zhou, H. (2021) Thermal Effect of a Revolving Gaussian Beam on Activating Heat-Sensitive Nociceptors in Skin. Journal of Applied Mathematics and Physics, 9. https://doi.org/10.4236/jamp.2021.91007
[8]	Cazares, S.M., Snyder, J.A., Belanich, J., Biddle, J.C., Buytendyk, A.M., Teng, S.H.M. and O’Connor, K. (2019) Active Denial Technology Computational Human Effects End-to-End Hypermodel for Effectiveness (ADT CHEETEH-E). Human Factors and Mechanical Engineering for Defense and Safety, 3, Article No. 13. https://doi.org/10.1007/s41314-019-0023-7
[9]	Parker, J.E., Nelson, E.J., Beason, C.W. and Cook, M.C. (2017) Effects of Variable Spot Size on Human Exposure to 95-GHz Millimeter Wave Energy. Technical Report, AFRL-RH-FS-TR-2017-0017.