Minimum Energy Requirement for Inducing Withdrawal Reflex in Millimeter Wave Exposures

Hongyun Wang^1*, Wesley A. Burgei², Shannon E. Foley², Hong Zhou³

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

DOI: 10.4236/jamp.2022.107162                

Abstract

We consider the problem of inducing withdrawal reflex on a test subject via exposure to a millimeter wave beam. In our physical model, there are 10 physical parameters affecting the occurrence of withdrawal reflex. Our goal is to pinpoint the roles of these physical parameters in inducing withdrawal reflex. We first carry out non-dimensionalization to reduce the model to a non-dimensional system of only 3 composite parameters: non-dimensional beam power density, non-dimensional beam radius, and non-dimensional exposure time. If the beam power is kept on and steady, withdrawal reflex occurs eventually; the shortest exposure time for inducing withdrawal reflex corresponds to the smallest energy consumption at the given power density and beam radius. In the 2D space of power density and beam radius, the overall minimum energy occurs at the corner of very large power density and very small beam radius, which also produces a very large value of maximum skin temperature and a long time to withdrawal reflex. To reduce the burn injury risk, we introduce a cap on the maximum skin temperature. At each given total beam power, we carry out optimizations with respect to the beam radius, constrained by the prescribed temperature cap. The energy consumption varies negatively with the prescribed temperature cap: a lower temperature cap can be accommodated only with a higher energy consumption via the venue of a larger beam radius. The energy consumption is relatively flat with respect to the total beam power and attains a minimum at a moderately large total beam power. The time to withdrawal reflex is approximately inversely proportional to the total beam power. Our analysis demonstrates that a moderately large total beam power is a good compromise to achieve both low energy consumption and short time to withdrawal reflex.

Keywords

Thermal Effect of Millimeter Wave in Skin, Activation of Thermal Nociceptors, Non-Dimensionalization Analysis, Optimization

1. Introduction

Millimeter waves (MMW) are radio waves in the frequencies range of 30 - 300 gigahertz (GHz). When human body is exposed to MMW radiation, a fraction of the electromagnetic energy is absorbed into the skin and it increases the skin temperature [1] [2]. The thermal nociceptors are activated once the activation temperature is reached, which leads to a heating sensation. As more nociceptors are activated, the heating sensation grows stronger and eventually it induces withdrawal reflex on the test subject and compels them to move away from the MMW beam [1] [2].

We consider the problem of inducing withdrawal reflex on a test subject while 1) minimizing the energy consumption and 2) controlling the maximum skin temperature to reduce the burn injury risk. For the beam energy absorption and the skin temperature evolution, we assume a semi-infinite skin tissue of uniform material properties (a single skin layer) and a uniform initial (baseline) skin temperature [3] [4]. We further assume that nociceptors are uniformly distributed in skin so that the heating sensation is quantitatively described by the skin volume in which nociceptors are activated (the activated volume). The occurrence of withdrawal reflex is determined by comparing the current activated volume and a threshold value on the volume. These model components, when assembled together, form a system that mathematically predicts the exposure time needed for inducing withdrawal reflex, given the MMW beam specifications [3] [4]. As a practical measure of energy consumption, we examine the electromagnetic energy absorbed into the skin, which has the advantage of being readily available in the model system above. To facilitate the theoretical analysis and numerical simulations, the model system is non-dimensionalized, which drastically reduces it to a system of 3 composite parameters: the non-dimensional beam power density, the non-dimensional beam radius, and the non-dimensional exposure time. In the non-dimensional system, the occurrence of withdrawal reflex is completely determined by these 3 composite parameters. In addition, the activated volume is proportional to the beam radius squared, which reduces its computational complexity to a function of 2 parameters.

In the non-dimensional system, we carry out two separate constrained minimizations of energy consumption with respect to the beam radius, subject to a prescribed maximum skin temperature. The two constrained optimizations are respectively, 1) at each given value of beam center power density, 2) at each given value of total beam power.

The rest of the paper is organized as follows. In Section 2, we describe the assumptions and set up the mathematical formulations of the model. Non-dimensionalization is carried out in Section 3, to reduce the number of parameters in the system. In Section 4, we explore the phenomenon and mechanism of continued rising of activated volume beyond the exposure time (i.e., after the beam power is turned off). With continued rising of activated volume, withdrawal reflex may occur at a time after the beam power is turned off. The minimum exposure time for inducing withdrawal reflex, at each given pair of power density and beam radius, is studied in Section 5. For a given beam, the minimum exposure time corresponds to the smallest energy consumption. In Section 6, at either a fixed beam power density or a fixed total beam power, we tune the beam radius to minimize the energy consumption while complying with the prescribed cap on skin maximum temperature. In Section 7, useful empirical formulas are constructed for the conditional minimization of energy consumption. We discuss and summarize the main findings in Section 8.

2. Mathematical Model

We consider the situation where a skin area of a test subject is exposed to an electromagnetic beam of millimeter wavelength [5] [6] [7]. We adopt a model similar to the one in our previous studies [3] [4] [8], which we now describe briefly. The model is based on the assumptions below.

1) The skin surface exposed to the millimeter wave is flat, which we select as the xy plane of the coordinate system.

2) The beam ray is perpendicular to the skin surface and is selected as the z axis.

3) Over the xy plane, the power density has an axisymmetric normal distribution around the beam center, which we set as the origin of the coordinate system.

$P\left(r\right)=P_d^{\left(i\right)}\exp \left(-2\frac{{\left|r\right|}^2}{r_{\text{b}}^2}\right)\equiv P_d^{\left(i\right)}G\left(r\right)\mathrm{,}r\equiv \left(x\mathrm{,}y\right)$ (1)

In (1), $P_d^{\left(i\right)}$ is the incident beam power density at the center and $r_{\text{b}}$ is the radius of the Gaussian beam [9]. $r_{\text{b}}$ is related to the radius of half-maximum region $r_{\text{HM}}$ by

$2\frac{r_{\text{HM}}^2}{r_{\text{b}}^2}=\ln 2$

$r_{\text{b}}=\sqrt{\frac{2}{\ln 2}}{r}_{\text{HM}}\approx 1.6986r_{\text{HM}}$ (2)

4) The skin’s material properties are uniform in $\left(z\mathrm{,}r\right)$.

5) Only the heat conduction along the z-direction (the depth direction) is included in the model; the lateral heat conduction in $r$ is neglected. This assumption is justified since the characteristic spatial scales in the lateral directions are much larger than the penetration depth of millimeter wave into the skin (less than 0.5 mm) [10].

6) The cooling effect of blood flow is neglected. The convection and radiation heat loss at the skin surface is also neglected. These simplifications are justified when the beam power is high and the exposure time is short [1] [2].

7) Before the start of exposure, the skin’s initial temperature is uniform in $\left(z\mathrm{,}r\right)$, which is called the skin baseline temperature and is denoted by $T_{\text{base}}$.

8) The activation of thermal nociceptors is governed by the activation temperature $T_{\text{act}}$. If a local temperature is above $T_{\text{act}}$, the nociceptors in that local region are activated.

9) The occurrence of withdrawal reflex is determined by the total number of nociceptors activated [11]. We consider the simple case where the spatial density of nociceptors is uniform in $\left(z\mathrm{,}r\right)$. The number of activated nociceptors is proportional to the activated volume. When the activated volume reaches a threshold value $v_c$, withdrawal reflex occurs; the test subject turns off the beam power and/or moves away to avoid the beam radiation [1] [2].

Let $t=0$ be the beam start time, and $T\left(r\mathrm{,}z\mathrm{,}t\right)$ the skin temperature at lateral location $r$, depth z, and time t. The temperature distribution $T\left(r\mathrm{,}z\mathrm{,}t\right)$ is governed by

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

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(r\right)$ is the Gaussian distribution of beam power density in (1),

· $\alpha$ is the fraction of beam power absorbed into the skin, and

· $t_{\text{end}}$ is the beam end time ( $t=0$ is set to the beam start time).

In (3), $1_{\left[\mathrm{0,}{t}_{\text{end}}\right]}\left(t\right)$ is the indicator function for the exposure period $\left[\mathrm{0,}{t}_{\text{end}}\right]$.

$1_{\left[\mathrm{0,}{t}_{\text{end}}\right]}\left(t\right)=\{\begin{array}{ll}\mathrm{1,}\hfill & t\in \left[\mathrm{0,}{t}_{\text{end}}\right]\hfill \\ \mathrm{0,}\hfill & \text{otherwise}\hfill \end{array}$ (4)

Initial boundary value problem (3) does not involve any derivative with respect to the lateral location $r$. The only influence of $r$ is in the coefficient $P\left(r\right)$ of the heat source term. As a result, we can separate out the dependence on $r$. We write $T\left(r\mathrm{,}z\mathrm{,}t\right)$ as

$T\left(r\mathrm{,}z\mathrm{,}t\right)=T_{\text{base}}+\left(\alpha P_d^{\left(i\right)}\right)G\left(r\right)U\left(z\mathrm{,}t\right)$

where $U\left(z\mathrm{,}t\right)$ is the temperature increase per absorbed power density and is governed by

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

The activated volume at time t is

$\begin{array}{c}{V}_{\text{act}}\left(t\right)=\text{Volume}\left\{\left(r\mathrm{,}z\right)|T\left(r\mathrm{,}z\mathrm{,}t\right)\ge T_{\text{act}}\right\}\\ =\text{Volume}\left\{\left(r\mathrm{,}z\right)|\left(\alpha P_d^{\left(i\right)}\right)G\left(r\right)U\left(z\mathrm{,}t\right)\ge \left(T_{\text{act}}-T_{\text{base}}\right)\right\}\end{array}$ (6)

Withdrawal reflex occurs when the activated volume reaches the threshold value, $v_c$. For an exposure event with a prescribed exposure time, it results in withdrawal reflex if ${\max }_t{V}_{\text{act}}\left(t\right)\ge v_c$. As we will see, ${\max }_t{V}_{\text{act}}\left(t\right)$ may be attained at a time $t>t_{\text{end}}$ (i.e., after the beam power is turned off).

In the definition of activated volume in (6), on the left side of the inequality condition, the multiplier ( $\alpha P_d^{\left(i\right)}$ ) is the beam power density absorbed into the skin, the distribution of relative beam power density $G\left(r\right)$ contains the effect of beam radius $r_{\text{b}}$, and solution $U\left(z\mathrm{,}t\right)$ of (5) contains the effects of exposure time $t_{\text{end}}$ and skin material properties $\left({\rho }_{\text{m}}\mathrm{,}{C}_p\mathrm{,}k\mathrm{,}\mu \right)$. The right side of the inequality condition contains the skin baseline temperature $T_{\text{base}}$ and the nociceptor activation temperature $T_{\text{act}}$. Withdrawal reflex occurs when the activated volume reaches the threshold value $v_c$. In total, there are 10 physical quantities affecting the occurrence of withdrawal reflex.

$\left(\alpha P_d^{\left(i\right)}\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)\mathrm{,}\left({\rho }_{\text{m}}\mathrm{,}{C}_p\mathrm{,}k\mathrm{,}\mu \right)\mathrm{,}\left(T_{\text{base}}\mathrm{,}{T}_{\text{act}}\right)\mathrm{,}{v}_c$ (7)

To simplify the mathematical formulation and to pinpoint the roles of these physical quantities in inducing withdrawal reflex, we carry out non-dimensionalization.

3. Non-Dimensionalization and Analytical Solution

For each physical quantity, we introduce a scale for non-dimensionalization. The depth scale is provided by, 1/μ, the characteristic penetration depth of millimeter wave into the skin. The length scale in the lateral directions is derived from volume threshold $v_{\text{c}}$.

· Length scale in the depth direction and non-dimensional depth

$z_{\text{s}}=\frac{1}{\mu }\mathrm{,}{z}_{\text{nd}}=\frac{z}{z_{\text{s}}}$

· Length scale in the lateral directions and non-dimensional lateral coordinates

$r_{\text{s}}=\sqrt{\frac{v_{\text{c}}}{\pi z_{\text{s}}}}=\sqrt{\frac{\mu v_{\text{c}}}{\pi }}\mathrm{,}{r}_{\text{nd}}=\frac{r}{r_{\text{s}}}$

Here $r_{\text{s}}$ is defined as the radius of the cylinder with volume $v_c$ and height $z_{\text{s}}$.

· Volume scale and non-dimensional volume

$v_{\text{s}}=r_{\text{s}}^2{z}_{\text{s}}=\frac{v_{\text{c}}}{\pi }\mathrm{,}{V}_{\text{nd}}=\frac{V}{v_{\text{s}}}$

· Time scale and non-dimensional time

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

· Temperature scale and non-dimensional temperature

$\Delta T_{\text{s}}=\left(T_{\text{act}}-T_{\text{base}}\right)\mathrm{,}{T}_{\text{nd}}=\frac{T-T_{\text{base}}}{\Delta T_{\text{s}}}$

· Power density scale and non-dimensional power density

$P_{\text{s}}=\frac{z_s{\rho }_{\text{m}}{C}_p\Delta T_{\text{s}}}{t_s}=k\mu \Delta T_{\text{s}}\mathrm{,}{P}_{\text{nd}}=\frac{P}{P_{\text{s}}}$

· Scale for $U\equiv \left(T-T_{\text{base}}\right)/\left(\alpha P\right)$ and non-dimensional U

$U_{\text{s}}=\frac{\Delta T_{\text{s}}}{P_{\text{s}}}=\frac{1}{k\mu }\mathrm{,}{U}_{\text{nd}}=\frac{U}{U_{\text{s}}}$

· Energy scale and non-dimensional energy

$E_{\text{s}}=P_{\text{s}}{r}_{\text{s}}^2{t}_{\text{s}}=\frac{{\rho }_{\text{m}}{C}_p{v}_c}{\pi }\Delta T_{\text{s}}\mathrm{,}{E}_{\text{nd}}=\frac{E}{E_{\text{s}}}$

For simplicity, we shall drop the subscript nd in all non-dimensional quantities. Instead, the original physical quantities will be distinguished with subscript phy when necessary for clarity. For example, $P_d^{\left(i\right)}$ refers to the non-dimensional incident beam power density at beam center while $P_{\text{d,phy}}^{\left(i\right)}$ is the physical incident beam power density before non-dimensionalization. With this notation, the non-dimensional versions of baseline temperature, activation temperature and volume threshold are

$T_{\text{base}}=0,T_{\text{act}}=\mathrm{1,}{v}_c=\pi$ (8)

The non-dimensional beam power density distribution is $\left(\alpha P_d^{\left(i\right)}\right)G\left(r\mathrm{<mo>\; ;\; </mo>}{r}_{\text{b}}\right)$ where

$r=r_{\text{phy}}\sqrt{\frac{\pi }{\mu v_{\text{c,phy}}}}\mathrm{,}{r}_{\text{b}}=r_{\text{b,phy}}\sqrt{\frac{\pi }{\mu v_{\text{c,phy}}}}$ (9)

$P_d^{\left(i\right)}=\frac{P_{\text{d,phy}}^{\left(i\right)}}{k\mu \left(T_{\text{act,phy}}-T_{\text{base,phy}}\right)}$ (10)

$G\left(r\mathrm{<mo>\; ;\; </mo>}{r}_{\text{b}}\right)=\exp \left(-2\frac{{\left|r\right|}^2}{r_{\text{b}}^2}\right)$

The non-dimensional version of $U\left(z\mathrm{,}t\right)$ is governed by

$\{\begin{array}{l}\frac{\partial U}{\partial t}=\frac{{\partial }^{2}U}{\partial z^2}+1_{\left[\mathrm{0,}{t}_{\text{end}}\right]}\left(t\right)\exp \left(-z\right)\hfill \\ {\frac{\partial U}{\partial z}|}_{z\mathrm{<mo>\; =\; </mo>0}}=\mathrm{0,}{U|}_{t\mathrm{<mo>\; =\; </mo>0}}=0\hfill \end{array}$ (11)

Notice that function $U\left(z\mathrm{,}t\right)$ has only one parameter, the non-dimensional beam end time

$t_{\text{end}}=t_{\text{end,phy}}\frac{k{\mu }^2}{{\rho }_{\text{m}}{C}_p}$ (12)

The non-dimensional temperature distribution is

$T\left(r\mathrm{,}z\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}\left(\alpha P_d^{\left(i\right)}\right)\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)=\left(\alpha P_d^{\left(i\right)}\right)G\left(r\mathrm{<mo>\; ;\; </mo>}{r}_{\text{b}}\right)U\left(z\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{t}_{\text{end}}\right)$

Absorption fraction $\alpha$ and incident beam power density $P_d^{\left(i\right)}$ appear only in the combination ( $\alpha P_d^{\left(i\right)}$ ). For conciseness, we introduce the power density absorbed into the skin:

$P_d\equiv \left(\alpha P_d^{\left(i\right)}\right)=\frac{\alpha P_{\text{d,phy}}^{\left(i\right)}}{k\mu \left(T_{\text{act,phy}}-T_{\text{base,phy}}\right)}$ (13)

An exposure event is completely specified by parameters $\left(P_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)$. The corresponding non-dimensional activated volume at time t is

$V_{\text{act}}\left(t\mathrm{<mo>\; ;\; </mo>}{P}_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)=\text{Volume}\left\{\left(r\mathrm{,}z\right)|P_{d}G\left(r\mathrm{<mo>\; ;\; </mo>}{r}_{\text{b}}\right)U\left(z\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{t}_{\text{end}}\right)\ge 1\right\}$ (14)

A change of variables $r_{\text{new}}=r_{\text{old}}/r_{\text{b}}$ separates out the dependence on $r_{\text{b}}$.

$\begin{array}{c}{V}_{\text{act}}\left(t\mathrm{<mo>\; ;\; </mo>}{P}_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)=r_{\text{b}}^2\text{Volume}\left\{\left(r\mathrm{,}z\right)|{P_{d}G\left(r\right)|}_{r_{\text{b}}=1}U\left(z\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{t}_{\text{end}}\right)\ge 1\right\}\\ =r_{\text{b}}^2{V_{\text{act}}\left(t\mathrm{<mo>\; ;\; </mo>}{P}_d\mathrm{,}{t}_{\text{end}}\right)|}_{r_{\text{b}}=1}\end{array}$ (15)

For a prescribed set of $\left(P_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)$, the exposure event leads to withdrawal reflex if

$\underset{t}{\max }{V}_{\text{act}}\left(t\mathrm{<mo>\; ;\; </mo>}{P}_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)\ge \pi$ (16)

where the maximum may be attained at $t>t_{\text{end}}$. Condition (16) is completely specified by 3 composite parameters $\left(P_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)$, defined in (9), (12) and (13).

The scaling property of the activated volume with respect to the beam radius in (15) is important for efficient numerical computation. It reduces the computational complexity of $V_{\text{act}}\left(t\mathrm{<mo>\; ;\; </mo>}{P}_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)$ to that of ${V_{\text{act}}\left(t\mathrm{<mo>\; ;\; </mo>}{P}_d\mathrm{,}{t}_{\text{end}}\right)|}_{r_{\text{b}}=1}$, a 2-parameter function.

One advantage of using the single layer skin model of uniform material properties in formulation (11) is that the solution $U\left(z\mathrm{,}t\right)$ has a closed-form analytical expression [12].

$U\left(z\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{t}_{\text{end}}\right)=u\left(z\mathrm{,}t\right)-u\left(z\mathrm{,}t-t_{\text{end}}\right)$ (17)

$u\left(z\mathrm{,}t\right)=\{\begin{array}{ll}-{\text{e}}^{-z}-z\text{erfc}\frac{z}{2\sqrt{t}}+\frac{2\sqrt{t}}{\sqrt{\pi }}{\text{e}}^{\frac{-z^2}{4t}}+\frac{{\text{e}}^{t-z}}{2}\text{erfc}\frac{2t-z}{2\sqrt{t}}+\frac{{\text{e}}^{t+z}}{2}\text{erfc}\frac{2t+z}{2\sqrt{t}}\mathrm{,}\hfill & t>0\hfill \\ \mathrm{0,}\hfill & t\le 0\hfill \end{array}$

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$ (18)

The time derivative of $u\left(z\mathrm{,}t\right)$ satisfies

$\frac{\partial u\left(z\mathrm{,}t\right)}{\partial t}=\frac{{\text{e}}^{t-z}}{2}\text{erfc}\frac{2t-z}{2\sqrt{t}}+\frac{{\text{e}}^{t+z}}{2}\text{erfc}\frac{2t+z}{2\sqrt{t}}>0$ (19)

It follows that for $t<t_{\text{end}}$ (before the beam power ends), $U\left(z\mathrm{,}t\right)=u\left(z\mathrm{,}t\right)$ is an increasing function of t at any z. But for $t>t_{\text{end}}$ (after the beam power ends), $U\left(z\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{t}_{\text{end}}\right)=u\left(z\mathrm{,}t\right)-u\left(z\mathrm{,}t-t_{\text{end}}\right)$ may decrease or increase with t, depending on the depth location z. In the next section, we study the temporal trend of $U\left(z\mathrm{,}t\mathrm{<mo>\; ;\; </mo>}{t}_{\text{end}}\right)$ for $t>t_{\text{end}}$, which is affected by $\left(z\mathrm{,}{t}_{\text{end}}\right)$, and we study the temporal trend of activated volume $V_{\text{act}}\left(t\right)$ which is affected by $\left(P_d\mathrm{,}{t}_{\text{end}}\right)$.

4. Continued Rising of Activated Volume after Beam Power Is Turned Off

In scaling law (15), the activated volume is proportional to the beam radius squared ( $r_{\text{b}}^2$ ). For the purpose of examining the temporal trend of activated volume, we only need to study the case of $r_b=1$, that is, we study ${V_{\text{act}}\left(t\mathrm{<mo>\; ;\; </mo>}{P}_d\mathrm{,}{t}_{\text{end}}\right)|}_{r_{\text{b}}=1}$ as a function of t.

At any given $\left(P_d\mathrm{,}{t}_{\text{end}}\right)$, the activated volume is initially zero for small t when the skin temperature is everywhere below the activation temperature. When the spatial maximum temperature (occurring at the beam center on the skin surface) reaches the activation temperature, $V_{\text{act}}\left(t\right)$ starts to increase from 0. While the beam power is kept on and steady, the activated $V_{\text{act}}\left(t\right)$ increases monotonically with t since the temperature $U\left(z\mathrm{,}t\right)$ increases with t at every depth location, as we discussed in (19). When the beam power is turned off, the temporal evolution of $U\left(z\mathrm{,}t\right)$ in (11) is driven solely by heat conduction. Heat flows from the high temperature region near the skin surface to the low temperature region inside skin. Figure 1 compares the non-dimensional temperature per power density $U\left(z\mathrm{,}t\right)$ vs z at $t=t_{\text{end}}$ and at $t=1.2t_{\text{end}}$ (a short time after the end of beam power). The heat flow at depth z is proportional to the spatial derivative $\partial U\left(z\mathrm{,}t\right)/\partial z$, which is zero at $z=0$, increases in magnitude with z for small z, and then decreases in magnitude with z for large z (see Figure 1). As a result, upon the beam power is turned off, for an interval at small depth z, there is a net heat out-flow and the temperature decreases with time t; for an interval at large depth z, there is a net heat in-flow and the temperature keeps increasing with time t even though the active heat source is off. In Figure 1, there are two regions in the z-direction, separated by $z_c\left(t_{\text{end}}\right)$, with opposite behaviors:

· for small depth $z<z_c\left(t_{\text{end}}\right)$, we have $U\left(z\mathrm{,1.2}{t}_{\text{end}}\right)<U\left(z\mathrm{,}{t}_{\text{end}}\right)$ ;

· for large depth $z>z_c\left(t_{\text{end}}\right)$, we have $U\left(z\mathrm{,1.2}{t}_{\text{end}}\right)>U\left(z\mathrm{,}{t}_{\text{end}}\right)$.

Here $z_c\left(t_{\text{end}}\right)$ is the location in z where ${\frac{\partial U\left(z\mathrm{,}t\right)}{\partial t}|}_{t_{\text{end}+\text{0}}}$ transitions from

negative to positive. The dependence of $z_c\left(t_{\text{end}}\right)$ on $t_{\text{end}}$ will be examined in Figure 3. We are interesting in the temporal trend of activated depth upon beam power is turned off. By definition, $d_{\text{act}}\left(t\right)$ satisfies $P_{d}U\left(d_{\text{act}}\left(t\right)\mathrm{,}t\right)=T_{\text{act}}=1$,

Figure 1. Evolution of $U\left(z\mathrm{,}t\right)$ upon beam power is turned off. For small z, $U\left(z\mathrm{,}t\right)$ decreases with t; for large z, $U\left(z\mathrm{,}t\right)$ increases with t. Accordingly, for low power density $P_d$, the activated depth $d_{\text{act}}\left(t\right)$ decreases with t; for high $P_d$, $d_{\text{act}}\left(t\right)$ increases with t.

which becomes

$U\left(d_{\text{act}}\left(t\right)\mathrm{,}t\right)=\frac{1}{P_d}$ (20)

In (20), the interaction between $1/P_d$ and the decreasing/increasing regions of $U\left(z\mathrm{,}t\right)$ vs t leads to decreasing/increasing trends of $d_{\text{act}}$ vs t, respectively, for low and for high $P_d$. In Figure 1, for a relatively low power density $P_{\text{d,1}}$, the value of $1/P_{\text{d,1}}$ is large, and the corresponding activated depth satisfies $d_{\text{act,1}}\left(t_{\text{end}}\right)<z_c\left(t_{\text{end}}\right)$. Recall that $U\left(z\mathrm{,1.2}{t}_{\text{end}}\right)<U\left(z\mathrm{,}{t}_{\text{end}}\right)$ for $z<z_c\left(t_{\text{end}}\right)$ and that $U\left(z\mathrm{,}t\right)$ is a decreasing function of z. To make $d_{\text{act,1}}\left(1.2t_{\text{end}}\right)$ satisfy (20) at $1.2t_{\text{end}}$, we must have $d_{\text{act,1}}\left(1.2t_{\text{end}}\right)<d_{\text{act,1}}\left(t_{\text{end}}\right)$. That is, for a relatively low power density $P_{\text{d,1}}$, upon the beam power is off, the activated depth starts decreasing. In contrast, for a relatively high power density $P_{\text{d,2}}$, the value of $1/P_{\text{d,2}}$ is small, and the corresponding activated depth satisfies $d_{\text{act,2}}\left(t_{\text{end}}\right)>z_c\left(t_{\text{end}}\right)$. Since $U\left(z\mathrm{,1.2}{t}_{\text{end}}\right)>U\left(z\mathrm{,}{t}_{\text{end}}\right)$ for $z>z_c\left(t_{\text{end}}\right)$ and $U\left(z\mathrm{,}t\right)$ decreases with z, (20) gives $d_{\text{act,2}}\left(1.2t_{\text{end}}\right)>d_{\text{act,2}}\left(t_{\text{end}}\right)$. That is, for a relatively high power density $P_{\text{d,2}}$, after the beam power is off, the activated depth will continue increasing for some time.

Similar to the situation with the activated depth $d_{\text{act}}\left(t\right)$ described above, the activated volume $V_{\text{act}}$ vs t also demonstrates opposite trends, respectively, for low $P_d$ and for high $P_d$. Let

$V_{\text{end}}\equiv V_{\text{act}}\left(t_{\text{end}}\right)\mathrm{,}{V}_{\text{max}}\equiv \underset{t}{\max }{V}_{\text{act}}\left(t\right)$ (21)

Figure 2 illustrates these two distinct trends of $V_{\text{act}}$ vs t. In the left panel,

Figure 2. Two distinct trends of $V_{\text{act}}$ vs t. Left panel: $P_d=2$ and $t_{\text{end}}=3$. Right panel: $P_d=8$ and $t_{\text{end}}=3$ ; $V_{\text{max}}$ is attained at $t_{\text{max}}>t_{\text{end}}$ with $V_{\text{max}}>V_{\text{end}}$.

$P_d=2$ and $t_{\text{end}}=3$ ; $V_{\text{act}}\left(t\right)$ decreases monotonically for $t>t_{\text{end}}$, and $V_{\text{end}}$ is the same as $V_{\text{max}}$. In the right panel, $P_d=8$ and $t_{\text{end}}=3$ ; $V_{\text{act}}\left(t\right)$ continues increasing for some time after the beam power is off; $V_{\text{max}}$ is attained at $t_{\text{max}}>t_{\text{end}}$ with $V_{\text{max}}>V_{\text{end}}$.

In Figure 1 and Figure 2, we see that at a fixed exposure time $t_{\text{end}}$ when the beam power density $P_d$ is large enough, the activated volume $V_{\text{act}}$ will continue rising after the beam power is off. In the discussion above we used the term “for a relatively large $P_d$ ” because in Figure 1 the continued rising of $V_{\text{act}}$ occurs when quantity ( $1/P_{\text{d}}$ ) is lower than the critical value $U\left(z_c\left(t_{\text{end}}\right)\mathrm{,}{t}_{\text{end}}\right)$ where $z_c\left(t_{\text{end}}\right)$ is the boundary in z between $U_t\left(z\mathrm{,}{t}_{\text{end}+\text{0}}\right)<0$ and $U_t\left(z\mathrm{,}{t}_{\text{end}+\text{0}}\right)>0$. In Figure 3, we show that both the boundary $z_c\left(t_{\text{end}}\right)$ and the critical value $U\left(z_c\left(t_{\text{end}}\right)\mathrm{,}{t}_{\text{end}}\right)$ are increasing functions of $t_{\text{end}}$. In particular, at any given beam power density $P_d$, when the exposure time $t_{\text{end}}$ is sufficiently large, $U\left(z_c\left(t_{\text{end}}\right)\mathrm{,}{t}_{\text{end}}\right)$ can be made larger than the given ( $1/P_d$ ) and consequently $V_{\text{max}}$ is attained at $t_{\text{max}}>t_{\text{end}}$ with $V_{\text{max}}>V_{\text{end}}$. In other words, any value of $P_d$ can qualify as a relatively high power density when $t_{\text{end}}$ is large enough. In summary, the continued rising of $V_{\text{act}}$ beyond $t_{\text{end}}$ can be realized at any given beam power density $P_d$ when the exposure time $t_{\text{end}}$ is large enough or at any given $t_{\text{end}}$ when $P_d$ is large enough.

The temporal trend of $V_{\text{act}}\left(t\right)$ is completely determined by two parameters $\left(P_d\mathrm{,}{t}_{\text{end}}\right)$. We examine the ratio $V_{\text{max}}/V_{\text{end}}$ as a function of $\left(P_d\mathrm{,}{t}_{\text{end}}\right)$. The top left panel of Figure 4 plots contours of $V_{\text{max}}/V_{\text{end}}$ vs $\left({\log }_2\left(P_d\right)\mathrm{,}{\log }_2\left(t_{\text{end}}\right)\right)$ ; the shaded region corresponds to $V_{\text{max}}=V_{\text{end}}$. The continued rising of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$, indicated by $V_{\text{max}}/V_{\text{end}}>1$, occurs when ${\log }_2\left(P_d\right)+{\log }_2\left(t_{\text{end}}\right)$ is large. The top right panel of Figure 4 shows contours of $d_{\text{act}}\left(t_{\text{end}}\right)$. As expected from the analysis in Figure 1, the continued rising of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$ occurs only when the activated depth is relatively large. This is most prominent at

Figure 3. $z_c\left(t_{\text{end}}\right)$ is the boundary between $U_t\left(z\mathrm{,}{t}_{\text{end}+\text{0}}\right)<0$ and $U_t\left(z\mathrm{,}{t}_{\text{end}+\text{0}}\right)>0$, introduced in Figure 1. Left: $z_c\left(t_{\text{end}}\right)$ vs $t_{\text{end}}$. Right: $U\left(z_c\left(t_{\text{end}}\right)\mathrm{,}{t}_{\text{end}}\right)$ vs $t_{\text{end}}$.

Figure 4. Comparison of contours of $V_{\text{max}}/V_{\text{end}}$ (top left), contours of $d_{\text{act}}\left(t_{\text{end}}\right)$ (top right), and contours of $T_{\text{max}}$ (bottom), in the 2D space of $\left({\log }_2\left(P_d\right)\mathrm{,}{\log }_2\left(t_{\text{end}}\right)\right)$.

low and moderate beam power densities. The bottom row of Figure 4 displays contours of the skin maximum temperature $T_{\text{max}}\equiv {\max }_{\left(r\mathrm{,}z\mathrm{,}t\right)}T\left(r\mathrm{,}z\mathrm{,}t\right)$. It is clear that the continued rising of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$ occurs only when the skin maximum temperature is very high. To get a 10% increase of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$ (i.e., $V_{\text{max}}/V_{\text{end}}=1.1$ ), the non-dimensional maximum temperature needs to reach $T_{\max }=8$ ; a 20% increase of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$ is associated with $T_{\text{max}}\approx 12$. Because of the very high skin temperature accompanying any meaningful increase of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$, it is impractical to realize and utilize such a gain.

There is another obstacle in harvesting the gain of continued rising of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$ for inducing withdrawal reflex: the optimal exposure design for such a purpose requires a fairly small beam radius. To harvest the gain and to maximize the energy efficiency, we need $V_{\text{max}}>V_{\text{end}}$ and we need to place the occurrence of withdrawal reflex right at $V_{\text{max}}$, which gives a constraint on the

optimal beam radius: ${r_{\text{b,gain}}^2{V}_{\text{max}}\left(P_d\mathrm{,}{t}_{\text{end}}\right)|}_{r_{\text{b}}=1}=\pi$. Here $r_{\text{b,gain}}$ is the optimal

beam radius for harvesting the gain of continued rising of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$ to induce withdrawal reflex. This constraint leads to

$r_{\text{b,gain}}\left(P_d\mathrm{,}{t}_{\text{end}}\right)={\left(\frac{\pi }{{V_{\text{max}}\left(P_d\mathrm{,}{t}_{\text{end}}\right)|}_{r_{\text{b}}=1}}\right)}^{1/2}$ (22)

$r_{\text{b,gain}}$ varies with $\left(P_d\mathrm{,}{t}_{\text{end}}\right)$, the prescribed power density and exposure time. Figure 5 shows contours of $r_{\text{b,gain}}$ in the 2D space of $\left({\log }_2\left(P_d\right)\mathrm{,}{\log }_2\left(t_{\text{end}}\right)\right)$ and compares them with those of $V_{\text{max}}/V_{\text{end}}$. The gains of 10% and 20% continued rising of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$ are associated with respectively $r_{\text{b,gain}}\le 0.8$ and $r_{\text{b,gain}}\le 0.65$ for the range of power density examined. For $r_{\text{b}}>1.5$, it is not possible to harvest any gain of continued rising of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$ for

Figure 5. Comparison of contours of $V_{\text{max}}/V_{\text{end}}$ and contours of $r_{\text{b,gain}}$ in the 2D space of $\left({\log }_2\left(P_d\right)\mathrm{,}{\log }_2\left(t_{\text{end}}\right)\right)$. Here $r_{\text{b,gain}}$ is the optimal beam radius for harvesting the gain in $V_{\text{max}}/V_{\text{end}}$ to induce withdrawal reflex.

inducing withdrawal reflex. Recall that continued rising of $V_{\text{act}}\left(t\right)$ beyond $t_{\text{end}}$ occurs when $t_{\text{end}}$ is sufficiently large. With a beam of radius $r_{\text{b}}>1.5$, withdrawal reflex would occur before the $t_{\text{end}}$ required for bringing on the continued rising of $V_{\text{act}}\left(t\right)$.

5. Energy Required for Inducing Withdrawal Reflex

At each given set of power density and beam radius $\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, when the exposure time $t_{\text{end}}$ is short enough, no volume is activated and no withdrawal reflex occurs. With the beam power density kept steady, as the exposure time increases, the temperature at each z increases unbounded, and thus, the activated volume increases unbounded. When the exposure time $t_{\text{end}}$ is long enough, we have ${\max }_t{V}_{\text{act}}\left(t\right)\ge \pi$ and the exposure event results in withdrawal reflex, which may occur after the beam power is turned off (i.e., $t_{\text{ref}}>t_{\text{end}}$ ). Let $t_{\text{end,min}}$ denote the minimum exposure time for inducing withdrawal reflex. At each given set of $\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, the total beam power is given and the minimum exposure time corresponds to the minimum energy consumption. With $t_{\text{end,min}}$, withdrawal reflex occurs at the temporal maximum of $V_{\text{act}}\left(t\right)$, which may occur after $t_{\text{end}}$. Let $t_{\text{ref,minE}}$ denote the time of withdrawal reflex with the minimum exposure time $t_{\text{end,min}}$. Both $t_{\text{end,min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$ and $t_{\text{ref,minE}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$ vary with $\left(P_d\mathrm{,}{r}_{\text{b}}\right)$. It should be pointed out that $t_{\text{ref,minE}}$ is the time of withdrawal reflex with the smallest energy consumption (i.e., the shortest exposure time), not necessarily the shortest time to withdrawal reflex. At a given set of $\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, if $t_{\text{ref,minE}}\mathrm{<mo>\; >\; </mo>}{t}_{\text{end,min}}$, then keeping the beam power on longer than $t_{\text{end,min}}$ will bring on withdrawal reflex sooner than $t_{\text{ref,minE}}$ at the expense of higher energy consumption.

To describe the mathematical formulation for $t_{\text{end,min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$ and $t_{\text{ref,minE}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, we introduce $V_{\text{max}}\left(P_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)$, the maximum activated volume achieved with the given $\left(P_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)$.

$V_{\text{max}}\left(P_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)\equiv \underset{t}{\max }{V}_{\text{act}}\left(t\mathrm{<mo>\; ;\; </mo>}{P}_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)=r_{\text{b}}^2\underset{t}{\max }{V_{\text{act}}\left(t\mathrm{<mo>\; ;\; </mo>}{P}_d\mathrm{,}{t}_{\text{end}}\right)|}_{r_{\text{b}}=1}$ (23)

From the definition, it is clear that $V_{\text{max}}\left(P_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)$ satisfies the scaling property

$V_{\text{max}}\left(P_d\mathrm{,}{r}_{\text{b}}\mathrm{,}{t}_{\text{end}}\right)=r_{\text{b}}^2{V_{\text{max}}\left(P_d\mathrm{,}{t}_{\text{end}}\right)|}_{r_{\text{b}}=1}$ (24)

$V_{\text{max}}$ may be attained at $t>t_{\text{end}}$ with $V_{\text{max}}>V_{\text{end}}$. Given $\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, the minimum exposure time $t_{\text{end,min}}$ and the corresponding reflex time $t_{\text{ref,minE}}$ have the expressions

$t_{\text{end,min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)=\min \left\{t_{\text{end}}|{V_{\text{max}}\left(P_d\mathrm{,}{t}_{\text{end}}\right)|}_{r_{\text{b}}=1}\ge \frac{\pi }{r_{\text{b}}^2}\right\}$ (25)

$t_{\text{ref,minE}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)=\text{arg}\underset{t}{\max }{V_{\text{act}}\left(t\mathrm{<mo>\; ;\; </mo>}{P}_d\mathrm{,}{t}_{\text{end,min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)\right)|}_{r_{\text{b}}=1}$ (26)

Given $\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, the minimum amount of absorbed energy for inducing withdrawal reflex is

$\begin{array}{c}{E}_{\text{min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)=P_d{r}_{\text{b}}^2\left({\displaystyle \int {G\left(r\right)|}_{r_{\text{b}}=1}\text{d}x\text{d}y}\right)t_{\text{end,min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)\\ =\frac{\pi }{2}{r}_{\text{b}}^2{P}_d{t}_{\text{end,min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)\end{array}$ (27)

Note that although $V_{\text{max}}$ is proportional to $r_{\text{b}}^2$, the effect of $r_{\text{b}}$ on $t_{\text{end,min}}$ in (25) is complicated because of the nonlinear dependence of $V_{\text{max}}$ on $t_{\text{end}}$. We use a 2D grid to represent the space of $\left({\log }_2\left(P_d\right)\mathrm{,}{\log }_2\left(r_{\text{b}}\right)\right)$. At each grid point, we calculate $t_{\text{end,min}}$, $E_{\text{min}}$, $t_{\text{ref,minE}}$, $T_{\text{max}}$ (the maximum skin temperature over space and over time), and $d_{\text{act}}\left(t_{\text{ref,minE}}\right)$ (the activated depth at reflex). Figure 6 compares the contours of $E_{\text{min}}$, $t_{\text{ref,minE}}$, $T_{\text{max}}$, and $d_{\text{act}}\left(t_{\text{ref,minE}}\right)$. In the top left panel, the energy consumption increases mainly with the beam radius. If energy efficiency is the only concern, low energy consumption may be achieved by using beams of small radius to induce withdrawal reflex. This approach, however, has several undesirable side effects

Figure 6. Comparison of contours of $E_{\text{min}}$ (top left), contours of $t_{\text{ref,minE}}$ (top right), contours of $T_{\text{max}}$ (bottom left), and contours of $d_{\text{act}}\left(t_{\text{ref,minE}}\right)$ (bottom right), in the 2D space of $\left({\log }_2\left(P_d\right)\mathrm{,}{\log }_2\left(r_{\text{b}}\right)\right)$.

· The time to withdrawal reflex is longer with smaller beam radius (top right panel of Figure 6). The longer time to reflex may not be desirable for the purpose of inducing withdrawal reflex as quickly as possible.

· When the beam radius is small, the maximum skin temperature is very large (bottom left panel of Figure 6). To prevent thermal injury, we need to set a cap on the maximum skin temperature allowed [13] [14] [15].

· Small beam radius leads to large activated depth at reflex (bottom right panel of Figure 6), which may not be suitable for the purpose of preventing thermal injury to the tissue deep inside skin.

Figure 7 displays $E_{\text{min}}$ and $T_{\text{max}}$ vs beam radius $r_{\text{b}}$ at several values of power density $P_d$. The energy required for inducing withdrawal reflex, $E_{\text{min}}$ is flat with respect to $r_{\text{b}}$ for small $r_{\text{b}}$ and starts increasing with $r_{\text{b}}$ for $r_{\text{b}}>1$. The maximum temperature $T_{\text{max}}$ is high for small $r_{\text{b}}$ and decreases rapidly toward 1 as $r_{\text{b}}$ is increased. Since $T_{\text{act}}=1$ (nondimensional), $T_{\max }=1$ is the lowest maximum skin temperature for having a positive activated volume and thus for inducing withdrawal reflex. Figure 8 examines $E_{\text{min}}$ and $T_{\text{max}}$ vs power density $P_d$ at several values of beam radius $r_{\text{b}}$. $E_{\text{min}}$ is relatively flat with respect to $P_d$ at each $r_{\text{b}}$ ; the height of flat $E_{\text{min}}$ vs $P_d$ increases significantly when $r_{\text{b}}$ is increased beyond $r_{\text{b}}=1$. For large $r_{\text{b}}$, $T_{\text{max}}$ is small (close to 1) and is relatively flat with respect to $P_d$. For small $r_{\text{b}}$, $T_{\text{max}}$ is large and increases with $P_d$.

6. Minimization of Energy Subject to a Prescribed Temperature Cap

From the top left panel of Figure 6, the unconstrained minimization of energy consumption with respect to $\left(P_d\mathrm{,}{r}_{\text{b}}\right)$ is attained at the smallest beam radius $r_{\text{b}}$ and the highest power density $P_d$ (the bottom right corner of the panel). The combination of small beam radius and large power density, however, is

Figure 7. $E_{\text{min}}$ and $T_{\text{max}}$ as functions of beam radius $r_{\text{b}}$ at several values of power density $P_d$. Both the horizontal and vertical axes have logarithmic scales.

Figure 8. $E_{\text{min}}$ and $T_{\text{max}}$ as functions of power density $P_d$ at several values of beam radius $r_{\text{b}}$. Both the horizontal and vertical axes have logarithmic scales.

accompanied by undesirable features of a very high maximum skin temperature $T_{\text{max}}$ (see the bottom left panel of Figure 6) and a very long time to reflex $t_{\text{ref,minE}}$ (see the top right panel of Figure 6). For practical applications, we consider the constrained optimization problem of minimizing energy consumption subject to $T_{\text{max}}\le T_{\text{cap}}$ where $T_{\text{cap}}$ is the prescribed cap on the maximum skin temperature.

Recall that at each given set of $\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, the shortest exposure time for inducing withdrawal reflex is $t_{\text{end,min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$ given in (25). Associated with this shortest exposure time, we have i) the energy consumption $E_{\text{min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, ii) the time of withdrawal reflect $t_{\text{ref,minE}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, which may be later than $t_{\text{end,min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, iii) the maximum skin temperature over space and time $T_{\text{max}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, and iv) the activated depth at withdrawal reflex $d_{\text{act}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$.

6.1. Minimization at a Given Beam Power Density P_d

At each given value of $P_d$, we optimize $E_{\text{min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)$ with respect to variable $r_{\text{b}}$, subject to the constraint $T_{\text{max}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)\le T_{\text{cap}}$. The conditionally optimal $r_{\text{b}}$ and the corresponding $E_{\text{min}}$, $t_{\text{end,min}}$, $t_{\text{ref,minE}}$, and $d_{\text{act}}$ are given by

$r_{\text{b}}\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]=\underset{r_{\text{b}}\mathrm{,}\text{s}\text{.t}\text{.}{T}_{\text{max}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)\le T_{\text{cap}}}{\text{argmin}}{E_{\text{min}}\left(P_d\mathrm{,}{r}_{\text{b}}\right)|}_{P_d\text{=given}}$ (28)

$E_{\text{min}}\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]=E_{\text{min}}\left(P_d\mathrm{,}{r}_{\text{b}}\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]\right)$

$t_{\text{end,min}}\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]=t_{\text{end,min}}\left(P_d\mathrm{,}{r}_{\text{b}}\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]\right)$

$t_{\text{ref,minE}}\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]=t_{\text{ref,minE}}\left(P_d\mathrm{,}{r}_{\text{b}}\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]\right)$

$d_{\text{act}}\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]=d_{\text{act}}\left(P_d\mathrm{,}{r}_{\text{b}}\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]\right)$

In Figure 9, we plot the optimal $r_{\text{b}}\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]$ and the corresponding $E_{\text{min}}$, $t_{\text{end,min}}$, $t_{\text{ref,minE}}$, and $d_{\text{act}}$ as functions of $P_d$ at $T_{\text{cap}}=\left\{\mathrm{1.5,2.0,2.5,3.0,3.5}\right\}$.

Figure 9. Conditional minimization of $E_{\text{min}}$ at each given power density $P_d$, subject to $T_{\text{max}}\le T_{\text{cap}}$. Plots of optimal $\left(r_{\text{b}}\mathrm{,}{E}_{\text{min}}\mathrm{,}{t}_{\text{end,min}}\mathrm{,}{t}_{\text{ref,minE}}\mathrm{,}{d}_{\text{act}}\right)$ vs $P_d$.

These are the results of conditional optimization given $\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]$. At any given $P_d$, when the prescribed cap on the maximum skin temperature $T_{\text{cap}}$ is lower, the optimal beam radius $r_{\text{b}}$ is larger (top left panel); the energy consumption $E_{\text{min}}$ is higher (top right panel); both the exposure time $t_{\text{end,min}}$ and the time to withdrawal reflect $t_{\text{ref,minE}}$ are lower (bottom left panel); and the activated depth at reflex $d_{\text{act}}$ is smaller (bottom right panel). In the bottom left panel, the exposure time $t_{\text{end,min}}$ is plotted as lines while the time to withdrawal reflect $t_{\text{ref,minE}}$ is represented by symbols. For $T_{\text{cap}}=3.5$ and $P_d>2^4=16$, we have $t_{\text{ref,minE}}>t_{\text{end,min}}$, implying that in this regime the gain of continued rising of $V_{\text{act}}$ beyond the beam end time is realized. Away from this regime, $t_{\text{ref,minE}}\mathrm{<mo>\; =\; </mo>}{t}_{\text{end,min}}$, indicating that it is impossible to realize the gain of continued rising of $V_{\text{act}}$ beyond the beam end time when both $P_d$ and $T_{\text{cap}}$ are moderate or small. At each $T_{\text{cap}}$, the conditionally optimal beam radius $r_{\text{b}}\left(P_d\right)$ is an increasing function of $P_d$ (top left panel); the corresponding energy consumption $E_{\text{min}}\left(P_d\right)$ exhibits a minimum with respect to $P_d$ at a moderate value of $P_d$ (top right panel); both the exposure time $t_{\text{end,min}}\left(P_d\right)$ and the time to withdrawal reflex $t_{\text{ref,minE}}\left(P_d\right)$ are approximately inversely proportional to $P_d$ (bottom left panel), and the activated depth at reflex $d_{\text{act}}\left(P_d\right)$ decreases with $P_d$ for small $P_d$ and is flat for large $P_d$ (bottom right panel). In summary, when the conditional minimization of energy is carried out at each given beam power density $P_d$, subject to constraint $T_{\text{max}}\le T_{\text{cap}}$, the minimum energy consumption $E_{\text{min}}$ is predominantly determined by the prescribed $T_{\text{cap}}$ while both the exposure time $t_{\text{end,min}}$ and the time to withdrawal reflect $t_{\text{ref,minE}}$ are mainly influenced by the beam power density $P_d$. Overall, a moderately large beam power density $P_d$ with the corresponding optimal beam radius satisfying the constraint $T_{\text{max}}\le T_{\text{cap}}$ is a good compromise to achieve 1) low energy consumption, 2) short exposure time and short time to withdrawal reflex, and 3) small activated depth at reflex.

6.2. Minimization at a Given Total Beam Power P_T

For a Gaussian beam with radius $r_{\text{b}}$, the beam center power density $P_d$ and the total beam power $P_{\text{T}}$ are related by

$P_{\text{T}}=\frac{\pi }{2}{r}_{\text{b}}^2{P}_d\mathrm{,}{P}_d=\frac{2P_{\text{T}}}{\pi r_{\text{b}}^2}$ (29)

At each given value of $P_{\text{T}}$, we optimize $E_{\text{min}}\left(\frac{2P_{\text{T}}}{\pi r_{\text{b}}^2}\mathrm{,}{r}_{\text{b}}\right)$ with respect to variable $r_{\text{b}}$, subject to the constraint $T_{\text{max}}\left(\frac{2P_{\text{T}}}{\pi r_{\text{b}}^2}\mathrm{,}{r}_{\text{b}}\right)\le T_{\text{cap}}$. The conditionally optimal

$r_{\text{b}}$ and the corresponding $P_d$, $E_{\text{min}}$, $t_{\text{end,min}}$, $t_{\text{ref,minE}}$, and $d_{\text{act}}$ are given by

$r_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=\underset{r_{\text{b}}\mathrm{,}\text{s}\text{.t}\text{.}{T}_{\text{max}}\left(\frac{2P_{\text{T}}}{\pi r_{\text{b}}^2}\mathrm{,}{r}_{\text{b}}\right)\le T_{\text{cap}}}{\text{argmin}}{E_{\text{min}}\left(\frac{2P_{\text{T}}}{\pi r_{\text{b}}^2}\mathrm{,}{r}_{\text{b}}\right)|}_{P_{\text{T}}\text{=given}}$ (30)

$P_d\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=\frac{2P_{\text{T}}}{\pi {\left(r_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}\right)}^2}$

$E_{\text{min}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=E_{\text{min}}\left(\frac{2P_{\text{T}}}{\pi {\left(r_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}\right)}^2}\mathrm{,}{r}_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}\right)$

$t_{\text{end,min}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=t_{\text{end,min}}\left(\frac{2P_{\text{T}}}{\pi {\left(r_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}\right)}^2}\mathrm{,}{r}_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}\right)$

$t_{\text{ref,minE}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=t_{\text{ref,minE}}\left(\frac{2P_{\text{T}}}{\pi {\left(r_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}\right)}^2}\mathrm{,}{r}_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}\right)$

$d_{\text{act}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=d_{\text{act}}\left(\frac{2P_{\text{T}}}{\pi {\left(r_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}\right)}^2}\mathrm{,}{r}_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}\right)$

For clarity, in (28) we used the notation $F\left[P_d\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right]$ to represent the results of minimization at a given beam power density $P_d$ whereas in (30) we used $F\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$ for the results of minimization at a given total beam power $P_{\text{T}}$.

In Figure 10, we plot the optimal $r_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$ and the corresponding $P_d$, $E_{\text{min}}$, $d_{\text{act}}$, $t_{\text{end,min}}$ and ( $t_{\text{ref,minE}}/t_{\text{end,min}}$ ) as functions of the total beam power $P_{\text{T}}$ at several values of prescribed temperature cap $T_{\text{cap}}=\left\{\mathrm{1.5,2.0,2.5,3.0,3.5}\right\}$. These are the results of conditional optimization given $\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$. The optimal beam radius $r_{\text{b}}$ increases with $P_{\text{T}}$ while decreasing as $T_{\text{cap}}$ is raised (top left panel). The corresponding beam power density $P_d$ increases both with $P_{\text{T}}$ and with $T_{\text{cap}}$ (top right panel). The energy consumption $E_{\text{min}}$ is relatively flat and exhibits a minimum with respect to $P_{\text{T}}$ while decreasing rapidly to a limit as $T_{\text{cap}}$ is raised (middle left panel). The activated depth at reflex $d_{\text{act}}$ decreases with $P_{\text{T}}$ while increasing with $T_{\text{cap}}$ (middle right panel). The minimum exposure time for inducing withdrawal reflex $t_{\text{end,min}}$ is approximately inversely proportional to $P_{\text{T}}$ while decreasing to a limit as $T_{\text{cap}}$ is raised (bottom left panel). In the bottom right panel, the ratio ( $t_{\text{ref,minE}}/t_{\text{end,min}}$ ) is above 1 only for $T_{\text{cap}}\ge 3.5$ and $P_{\text{T}}\ge 2^{4.6}=24.25$. This result indicates that the gain of continued rising of $V_{\text{act}}$ beyond the beam end time cannot be realized when both $T_{\text{cap}}$ and $P_{\text{T}}$ are moderate or small. Even when this gain is clearly realized at $T_{\text{cap}}=3.5$ and ${\log }_2\left(P_{\text{T}}\right)=7$, the amount of energy saved is fairly small (comparing $E_{\text{min}}$ for $T_{\text{cap}}=3.5$ and for $T_{\text{cap}}=3.0$ in the middle left panel). In summary, when the conditional minimization of energy is carried out at each given total beam power $P_{\text{T}}$, subject to constraint $T_{\text{max}}\le T_{\text{cap}}$, the minimum energy consumption $E_{\text{min}}$ is again predominantly determined by the prescribed cap on the maximum skin temperature $T_{\text{cap}}$. The optimal beam radius varies slowly with ${\log }_2\left(P_{\text{T}}\right)$ and as a result the corresponding power density $P_d$ is nearly proportional to $P_{\text{T}}$. Both the exposure time $t_{\text{end,min}}$ and the time to withdrawal reflect $t_{\text{ref,minE}}$ are nearly inversely proportional to the total beam power $P_{\text{T}}$. Overall, a moderately large or large $P_{\text{T}}$ with the corresponding optimal beam radius satisfying the constraint $T_{\text{max}}\le T_{\text{cap}}$ is a good compromise to achieve i) low energy consumption, ii) short exposure time and short time to withdrawal reflex, and iii) small activated depth at reflex.

7. Empirical Formulas for Optimal Quantities

The results presented in Sections 5 and 6 are obtained by analytically solving the skin temperature evolution model (11) and numerically computing the activated volume in (15). Based on the simulated results, we now construct empirical formulas for the conditional optimization at each given total power $P_{\text{T}}$, subject to the constraint $T_{\text{max}}\le T_{\text{cap}}$. Specifically, we first construct empirical formulas for the optimal beam radius $r_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$ and the minimum energy consumption

Figure 10. Conditional minimization of $E_{\text{min}}$ at each given total power $P_{\text{T}}$, subject to $T_{\text{max}}\le T_{\text{cap}}$. Plots of optimal $\left(r_{\text{b}}\mathrm{,}{P}_d\mathrm{,}{E}_{\text{min}}\mathrm{,}{d}_{\text{act}}\mathrm{,}{t}_{\text{end,min}}\mathrm{,}{t}_{\text{ref,minE}}/t_{\text{end,min}}\right)$ vs $P_{\text{T}}$.

$E_{\text{min}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$. Then we use these two formulas to derive the empirical formulas for the corresponding power density $P_d\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$ and the exposure time $t_{\text{end,min}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$.

In the dimension of $P_{\text{T}}$, we use a quadratic of ${\log }_2\left(P_{\text{T}}\right)$ as the function form for both $r_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$ and $E_{\text{min}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$.

$F\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=f_2\left(T_{\text{cap}}\right){\left({\log }_2{P}_{\text{T}}\right)}^2+f_1\left(T_{\text{cap}}\right)\left({\log }_2{P}_{\text{T}}\right)+f_0\left(T_{\text{cap}}\right)$ (31)

In the dimension of $T_{\text{cap}}$, the three coefficients in the quadratic form (31) are set to

$f_j\left(T_{\text{cap}}\right)=c_{j\mathrm{,0}}+\frac{c_{j\mathrm{,1}}}{{\left(T_{\text{cap}}-c_{j\mathrm{,2}}\right)}^{\alpha }}$ (32)

where $\alpha =1$ for $r_{\text{b}}$ and $\alpha =3/2$ for $E_{\text{min}}$. Fitting the simulated results of $r_{\text{b}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$ and $E_{\text{min}}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$ to the two-variable function form specified by (31) and (32), we obtain the empirical approximations below.

$\begin{array}{c}{r}_{\text{b}}^{\left(\text{emp}\right)}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=\left(-0.325+\frac{0.301}{T_{\text{cap}}-1}\right)\frac{{\left({\log }_2{P}_{\text{T}}\right)}^2}{100}\\ +\left(0.611+\frac{1.204}{T_{\text{cap}}-0.681}\right)\frac{{\log }_2{P}_{\text{T}}}{10}+\left(0.406+\frac{0.929}{T_{\text{cap}}-0.681}\right)\end{array}$ (33)

$\begin{array}{c}{E}_{\text{min}}^{\left(\text{emp}\right)}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=\left(1.305+\frac{15.19}{{\left(T_{\text{cap}}-0.848\right)}^{3/2}}\right)\frac{{\left({\log }_2{P}_{\text{T}}\right)}^2}{100}\\ +\left(-8.242+\frac{0.475}{{\left(T_{\text{cap}}-1\right)}^{3/2}}\right)\frac{{\log }_2{P}_{\text{T}}}{10}+\left(8.647+\frac{9.121}{{\left(T_{\text{cap}}-1\right)}^{3/2}}\right)\end{array}$ (34)

$P_d^{\left(\text{emp}\right)}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=\frac{P_{\text{T}}}{0.5\pi {\left(r_{\text{b}}^{\left(\text{emp}\right)}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}\right)}^2}$ (35)

$t_{\text{end,min}}^{\left(\text{emp}\right)}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}=\frac{E_{\text{min}}^{\left(\text{emp}\right)}\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}}{P_{\text{T}}}$ (36)

In Figure 11, we compare the optimal quantities $\left(r_{\text{b}}\mathrm{,}{E}_{\text{min}}\mathrm{,}{P}_d\mathrm{,}{t}_{\text{end,min}}\right)$ obtained in accurate numerical simulations and the corresponding quantities predicted using empirical formulas (33)-(36). These quantities are compared as functions of ${\log }_2\left(P_{\text{T}}\right)$ at several values of $T_{\text{cap}}$. Figure 11 demonstrates that the empirical approximations (33)-(36) (represented as symbols) match very well the accurate numerical results of the model system (plotted as lines) in the range of $\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$ considered.

To measure the accuracy of empirical formulas quantitatively, we examine the relative errors in these approximations. For quantity F, anyone of $\left(r_{\text{b}}\mathrm{,}{E}_{\text{min}}\mathrm{,}{P}_d\mathrm{,}{t}_{\text{end,min}}\right)$, we define the relative error of the empirical approximation $F^{\left(\text{emp}\right)}$ as

$\text{Err}\left\{F^{\left(\text{emp}\right)}\right\}\equiv \ln \left(F\right)-\ln \left(F^{\left(\text{emp}\right)}\right)\approx \frac{F-F^{\left(\text{emp}\right)}}{F}$ (37)

In Figure 12, we plot the relative errors of empirical approximations (33)-(36)

Figure 11. Comparison of the simulated $\left(r_{\text{b}}\mathrm{,}{E}_{\text{min}}\mathrm{,}{P}_d\mathrm{,}{t}_{\text{end,min}}\right)$ (plotted as lines) and the corresponding empirical approximations (33)-(36) (represented as symbols).

as functions of ${\log }_2\left(P_{\text{T}}\right)$ at several values of $T_{\text{cap}}$. The relative errors are well below 2% in the range of $\left\{P_{\text{T}}\mathrm{<mo>\; ;\; </mo>}{T}_{\text{cap}}\right\}$ considered.

8. Concluding Remarks

In this paper, we considered the problem of inducing withdrawal reflex on a test subject via exposure to a millimeter wave beam. If we keep the beam power on and steady, the volume of skin tissue where nociceptors are activated grows without bound, and eventually withdrawal reflex occurs when the heating sensation is above the subject’s tolerance.

The occurrence of withdrawal reflex is affected by 10 physical parameters of beam specifications, exposure time, skin material properties, skin baseline temperature, activation temperature of thermal nociceptors, and critical threshold on the activated volume for withdrawal reflex (see the list in (7)). To pinpoint the roles of these physical parameters, we carried out non-dimensionalization to reduce the model to a non-dimensional system of only 3 composite parameters:

Figure 12. Relative errors of empirical approximations (33)-(36), defined in (37), for optimal quantities $\left(r_{\text{b}}\mathrm{,}{E}_{\text{min}}\mathrm{,}{P}_d\mathrm{,}{t}_{\text{end,min}}\right)$.

non-dimensional (absorbed) beam power density, non-dimensional beam radius, and non-dimensional exposure time (see the model described in (14) and (16)). Another benefit of non-dimensionalization is that the non-dimensional temperature distribution is written in terms of a parameter-free function of $\left(z\mathrm{,}t\right)$ that has a simple closed-form expression (see the solution given in (17)).

With the 3-parameter non-dimensional system, we explored the phenomenon and mechanism of continued rising of activated volume beyond the beam end time. While the beam is on and steady, the temperature is monotonically increasing with time at all depth. After the beam power is turned off, near skin surface the temperature decreases with time while deep inside skin the temperature keeps increasing with time, driven by the heat flow from the hot region near skin surface. For a given beam end time, there is a critical depth separating the skin tissue into two regions: for depth smaller than the critical value, the temperature decreases with time after the beam power is off; for depth larger than the critical value, the temperature continues increasing with time after the beam end time (see Figure 1). If, by the beam end time, the activated depth has reached the critical depth, then the activated depth and the activated volume will continue increasing with time after the beam power is off. At any given beam end time, the activated depth increases with the beam power density. It follows that when the beam power density is sufficiently high, the activated volume will continue increasing with time after the beam end time (see Figure 2). We also demonstrated that at any given beam power density, with a sufficiently long exposure time, the activated volume will continue increasing with time after the beam end time (see the results in Figure 3). The continued rising of activated volume after the beam end time, however, is associated with very high values of maximum skin temperature (see Figure 4), which makes it impractical to utilize this gain due to the need of reducing the burn injury risk.

At each given set of $\left(P_d\mathrm{,}{r}_{\text{b}}\right)$, there is a minimum exposure time such that the maximum activated volume over time reaches the critical volume threshold for withdrawal reflex, which may occur after the beam end time. This minimum exposure time corresponds to the minimum energy consumption at the given $\left(P_d\mathrm{,}{r}_{\text{b}}\right)$. We found that the minimum energy consumption is relatively flat with respect to $P_d$ and increases significantly with $r_{\text{b}}$. Lower energy consumption is achieved with a smaller beam radius, which is associated with several undesirable features: high value of maximum skin temperature, long time to withdrawal reflex, large activated depth at reflex (see Figure 6 and Figure 7).

To find a compromise between reducing the energy consumption and preventing the maximum skin temperature from getting too high, we carried out two constrained minimizations of energy with respect to $r_{\text{b}}$, subject to

$T_{\text{max}}\le T_{\text{cap}}$, respectively, 1) at a given beam power density $P_d$ or 2) at a given total beam power $P_{\text{T}}$. The results of the two constrained minimizations depend on, respectively, $\left[P_d\mathrm{,}{T}_{\text{cap}}\right]$ or $\left\{P_{\text{T}}\mathrm{,}{T}_{\text{cap}}\right\}$. In the two constrained minimizations, the minimum energy consumption is relatively flat with respect to $P_d$ or $P_{\text{T}}$ ; it increases rapidly as the temperature cap $T_{\text{cap}}$ is lowered, indicating the trade-off that a lower value of maximum skin temperature can be achieved only with a larger beam radius at the price of a significantly higher energy consumption (see Figure 9 and Figure 10). When the beam power density $P_d$ is fixed, a lower value of $T_{\text{cap}}$ (achieved with a larger beam radius) yields a shorter time to withdrawal reflex (see Figure 9). In contrast, when the total beam power

$P_{\text{T}}=\frac{1}{2}\pi r_{\text{b}}^2{P}_d$ is fixed, a lower value of $T_{\text{cap}}$ (achieved with a larger beam radius)

yields a slightly longer time to withdrawal reflex (see Figure 10). When the temperature cap $T_{\text{cap}}$ is fixed, the time to withdrawal reflex is approximately inversely proportional to $P_d$ or $P_{\text{T}}$ in both constrained minimizations. The activated depth at reflex varies positively with the temperature cap $T_{\text{cap}}$ : a lower value of $T_{\text{cap}}$ yields a smaller activated depth at reflex. In addition, the activated depth at reflex varies negatively with $P_d$ or $P_{\text{T}}$ : at a fixed value of $T_{\text{cap}}$, the activated depth at reflex can be reduced significantly with a higher value of $P_d$ or $P_{\text{T}}$ (see Figure 9 and Figure 10). In the two constrained minimizations, the continued rising of activated volume beyond the beam end time does not occur at the optimum unless $T_{\text{cap}}$ is large or $P_d$ is large or both. In most of the parameter region examined, with the optimal beam radius, withdrawal reflex occurs before the onset of the continued rising of activated volume. As a result, in most of the parameter region examined in Figure 9 and Figure 10, withdrawal reflex occurs right at the beam end time (the minimum exposure time required for inducing withdrawal reflex). Given a prescribed cap on the maximum skin temperature $T_{\text{cap}}$, a moderately large or large total beam power with the corresponding optimal beam radius satisfying the constraint $T_{\text{max}}\le T_{\text{cap}}$ is a good compromise to achieve 1) low energy consumption, 2) short exposure time and short time to withdrawal reflex, and 3) small activated depth at reflex. The minimum energy consumption varies negatively with the prescribed $T_{\text{cap}}$. A lower value of $T_{\text{cap}}$ can be accommodated only with a higher energy consumption via the venue of a larger beam radius.

All conclusions obtained in this paper are based on accurately solving the skin temperature evolution model (11) and numerically computing the activated volume in (14). From the accurate numerical solutions, we constructed simple empirical formulas for the constrained optimization as functions of the given total beam power and the given temperature cap (see formulas (33)-(36)). The empirical formulas involve only rational functions of the input variables; they are easy to implement, require very little computing power, and are very accurate. These empirical formulas provide a quick and practical way of obtaining approximately the results of constrained optimizations without running full simulations.

Acknowledgements

The authors acknowledge the Joint Intermediate Force Capabilities Office of US 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 US Government.

Conflicts of Interest

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

References

[1]	Parker, J.E., Nelson, E.J. and Beason, C.W. (2017) Thermal and Behavioral Effects of Exposure to 30-kW, 95-GHz Millimeter Wave Energy. Technical Report, AFRL-RH-FS-TR-2017-0016.
[2]	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.
[3]	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
[4]	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
[5]	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
[6]	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
[7]	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
[8]	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, 88-100. https://doi.org/10.4236/jamp.2021.91007
[9]	Paschotta, R. (2008) “Effective Mode Area” in the Encyclopedia of Laser Physics and Technology. Wiley-VCH.
[10]	Durney, C.H., Massoudi, H. and Iskander, M.F. (1986) Radiofrequency Radiation Dosimetry Handbook. 4th Edition, USAF School of Aerospace Medicine, Brooks AFB.
[11]	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
[12]	Wang, H., Burgei, W.A. and Zhou, H. (2020) Analytical Solution of One-Dimensional Pennes’ Bioheat Equation. Open Physics, 18, 1084-1092. https://doi.org/10.1515/phys-2020-0197
[13]	Ozen, S., Helhel, S. and Bilgin, S. (2011) Temperature and Burn Injury Prediction of Human Skin Exposed to Microwaves: A Model Analysis. Radiation and Environmental Biophysics, 50, 483-489. https://doi.org/10.1007/s00411-011-0364-y
[14]	Ng, E.Y. and Chua, L.T. (2002) Prediction of Skin Burn Injury. Part 1: Numerical Modelling. Proceedings of the Institution of Mechanical Engineers, Part H, 216, 157-170. https://doi.org/10.1243/0954411021536379
[15]	Ng, E.Y. and Chua, L.T. (2002) Prediction of Skin Burn Injury. Part 2: Parametric and Sensitivity Analysis. Proceedings of the Institution of Mechanical Engineers, Part H, 216, 171-183. https://doi.org/10.1243/0954411021536388