Assessing skin thermal injury risk in tests of electromagnetic heating until flight

Abstract

We assess the risk of skin thermal injury in a scenario where a test subject is exposed to an electromagnetic beam until the occurrence of a flight action (repel). The physical process is modeled as follows. The absorbed electromagnetic power increases the skin temperature. When the temperature exceeds a certain threshold, thermal nociceptors are activated, generating a nerve signal. Once the activated skin volume reaches a threshold, the flight signal is triggered. After a delay corresponding to the human reaction time, the flight action is realized when the subject moves away. The injury risk is quantified using the thermal damage parameter calculated from the Arrhenius equation. This parameter depends on the power density absorbed into the skin, which is not directly measurable. In addition, the volume threshold for flight initiation is unknown. To circumvent these difficulties, we normalize the formulation and express the thermal damage parameter in terms of the occurrence time-of-flight action, which is reliably observed in exposure tests. This formulation of thermal injury provides a viable framework for investigating the effects of model parameters.

1. Introduction

Millimeter-wave (MMW) systems, operating in the 30–300 GHz range, are becoming increasingly important in both defense and civilian applications, such as wireless communications, security scanning, modern medicine, industrial heating, police traffic devices, and military technologies. These diverse applications prompt a better understanding of the biological effects of MMW exposures, which has been extensively studied by many authors. Some of these studies were compiled in a special issue of the Journal of Directed Energy.^1–8

One major concern with such systems is the safety of human exposure to MMWs. Earlier clinical and experimental studies on the biological effects of MMW exposure have found that the primary impact of MMWs is thermal, with the skin being the main organ absorbing the MMW energy.^9–11 The MMW energy penetrates less than 1 mm into the skin, causing the skin temperature to rise due to dielectric heating.

The local temperature rise of the skin depends on several factors, including the intensity and spot size of the MMW radiation, the duration of exposure, the absorption coefficient at the MMW frequency, the heat conductivity of the skin, heat exchange via blood perfusion, and heat loss at the skin surface (through black body radiation, air convection, sweating). ¹² Skin thermal damage occurs not only due to high temperature but also due to the duration of exposure to elevated temperatures. There are two potential thermal hazards from millimeter waves: thermal pain and burns. For example, in a 94-GHz MMW exposure with a 1.8-W/cm² beam, thermal pain begins when the surface temperature reaches approximately 44°C, which is about a 10°C increase over the typical baseline skin temperature.^13–15 In contrast, the threshold for burns corresponds to peak temperatures of approximately 60°C, requiring a beam energy fluence (energy delivered per unit area) that is two to three times higher than the fluence needed to produce thermal pain. ¹⁶

Most MMW systems operate at low power. However, certain commercial and military applications involve MMW systems that operate at high power. In these cases, there is a potential risk of accidental over-exposure to MMWs, which could lead to thermal injury. In the present study, we focus on the effects of a high-power MMW exposure on skin tissues in scenarios where the exposure continues until the occurrence of flight (repel). We investigate how factors such as human reaction time, MMW beam specifications, and other parameters influence skin thermal damage.

This paper is organized into four subsequent sections. Section 2 presents the mathematical formulation for modeling skin thermal damage caused by MMW radiation. Section 3 lists the numerical values of all parameters collected from published literature, except for the unknown volume threshold for flight initiation. To address this, we establish a normalized formulation that eliminates the need for the unknown volume threshold. Section 4 follows with numerical simulations, investigating the relationship between thermal damage and the observed time-of-flight action, as well as its sensitivity to model parameters. Finally, Section 5 discusses the results and concludes this paper.

2. Mathematical formulation

2.1. Overview of the model

We consider a scenario where the test subject’s skin is exposed to MMW radiation. The absorbed electromagnetic power acts as a heat source, raising the skin temperature. When the local temperature exceeds the activation threshold, thermal nociceptors are activated, transducing a nerve signal that is transmitted to the brain and contributes to the heat sensation. Once the heat sensation surpasses a subjective threshold, the brain sends a flight instruction signal via the nervous system to the muscles. Finally, the muscles initiate the flight response, causing the subject to move away from the electromagnetic beam.

We start by listing the assumptions used in formulating our model:

• The skin tissue is modeled as a single layer with uniform material properties, which are assumed to be independent of depth within the skin.

• All thermal nociceptors in the skin exhibit identical behavior, including (1) responding solely to the instantaneous temperature, (2) having the same deterministic activation temperature, and (3) upon activation by a local temperature exceeding the activation temperature, transducing an electrical signal of uniform magnitude.

• Thermal nociceptors are uniformly distributed within the skin, such that the total nociceptor signal is proportional to the activated skin volume.

• The heat sensation is proportional to the sum of nociceptor signals, which in turn is proportional to the activated skin volume. When the activated skin volume exceeds a subjective threshold, a flight signal is triggered.

• There is a time delay between the initiation of the flight signal and the execution of the flight action (i.e., when the subject begins to move). This delay is analogous to the simple reaction time (SRT) measured in psychophysical experiments.

• Once the subject begins to move, the beam radiation is considered terminated.

• The thermal damage is modeled using the Arrhenius equation.

The primary motivation for adopting these simplification assumptions is to build a rudimentary thermal injury model that is operational based on the currently available skin parameter values.

The flight is initiated when a sufficiently large nerve signal is transduced from the exposed skin spot. The flight action is eventually carried out after a time delay, which accounts for the nociceptor signal traveling from the exposed spot to the brain, the transduction of the flight instruction signal in the brain, and the flight signal traveling to the muscles. This time delay between flight initiation and flight action is referred to as the SRT or the human reaction time. The beam radiation on the skin is considered terminated only when the flight action occurs, meaning the test subject moves away from the beam. The beam radiation does not terminate at the point of flight initiation. Due to the delay in power-off, the human reaction time influences the thermal injury risk.

Flight initiation is determined by the aggregated nociceptor signal transduced from the exposed skin spot and a subjective threshold. In this study, we consider the idealized scenario where thermal nociceptors are uniformly distributed within the skin tissue. As a result, the aggregated nociceptor signal is proportional to the activated skin volume, where the local temperature exceeds the nociceptor activation threshold. Under this assumption, flight initiation is determined by the activated skin volume and its corresponding volume threshold.

During a high-power MMW exposure, while the beam power is on, the skin temperature increases monotonically due to the absorbed electromagnetic power, with the cooling effects of blood perfusion and surface heat loss being relatively negligible. As a result, the activated skin volume increases monotonically, and the flight initiation threshold is reached as the exposure continues. We study the risk of skin thermal damage in an experimental setting where the beam radiation continues until the occurrence of flight action. At the point of flight action, the active beam heating is turned off, and the spatial maximum skin temperature starts to decrease. However, the skin temperature remains elevated for some time after the beam power is terminated. To accurately assess the skin thermal damage risk, it is essential to monitor the temperature evolution well beyond the end of beam power, until the temperature returns to a level close to the baseline skin temperature. The skin temperature history is used to calculate the thermal damage parameter, Ω, based on an Arrhenius-type relation.^17,18

We need a measurable quantity to parameterize the thermal damage risk. The power density absorbed into the skin is a significant factor influencing the occurrence and severity of thermal injury. However, in exposure tests, the power density absorbed into the skin is not directly measured. In our current study, the process from temperature evolution to nociceptor activation, to flight initiation, and finally to flight action is modeled as a deterministic process. That is, there are no fluctuations in skin material properties, the thermal nociceptor activation temperature, the volume threshold for flight initiation, or the human reaction time. Under this model assumption, there is a one-to-one correspondence between the power density absorbed into the skin and the occurrence time-of-flight action. We propose two approaches.

1. We use the power density absorbed into the skin to parameterize the thermal damage risk. Although the absorbed power density is not directly measured in exposure tests, we can estimate it from the specified incident power density and an empirical absorption fraction coefficient.

2. We use the occurrence time-of-flight action to parameterize the thermal damage risk.

Under our current deterministic model, these two approaches are equivalent. Both approaches are affected if the deterministic model parameters are varied. Although the model is deterministic, we can apply it to different sets of parameter values to study the effect of variation in model parameters. This is done in the sensitivity study section of the current work. In this study, most results are presented using approach 2, while one figure uses approach 1. There are two reasons for this choice. (1) Reducing the occurrence time-of-flight action may be a goal in system design. Using the occurrence time-of-flight action to parameterize the injury risk directly informs the trade-off between reducing the occurrence time-of-flight action and minimizing the thermal injury risk. (2) The occurrence time-of-flight action is directly observed in exposure tests.

2.2. The coordinate system

We first establish the coordinate system as described in our previous study. ¹⁹ We consider the case of a flat skin surface, with the normal direction of the skin surface pointing into the skin tissue, selected as the positive $z$ -direction. The skin surface is defined at $z=0$ , and $z>0$ represents the skin tissue beneath the surface. In a plane perpendicular to the $z$ -axis (i.e., with $z=z_0$ ), the 2D coordinates are denoted by $\mathbf{r}=(x,y)$ . We consider the case where the incident electromagnetic beam is perpendicular to the skin surface (i.e., along the $z$ -direction), and the power density over its cross-section follows a Gaussian distribution, as given in Equation (9). The center of the incident beam corresponds to the location of the spatial maximum power density over the beam’s cross-section. On the skin surface, we set the beam center at the origin of the $\mathrm{xy}$ -plane, $\mathbf{r}=(0,0)$ . In the skin tissue, the 3D coordinates of a point are represented by $(\mathbf{r},z)=(x,y,z),z\ge 0$ .

2.3. Physical quantities and equations

We introduce the physical quantities and equations used in our model.

• $(\mathbf{r},z)=(x,y,z),z\ge 0$ denotes the 3D coordinates of a point in the skin tissue.

• $T(\mathbf{r},z,t)$ is the skin temperature at position $(\mathbf{r},z)$ at time $t$ .

• ${\rho }_m$ is the mass density, $C_p$ is the specific heat capacity, and $k$ is the heat conductivity of skin. In this study, we consider a single layer of uniform skin where all material properties are independent of position $(\mathbf{r},z)$ .

• $P_d\left(\mathbf{r}\right)$ is the power density passing through the skin surface at position $\mathbf{r}=(x,y)$ . In our model, $P_d\left(\mathbf{r}\right)$ is the power density absorbed into the skin tissue (which is converted to heat). It is given by the expression $P_d\left(\mathbf{r}\right)=\alpha P_d^{(i)}\left(\mathbf{r}\right)$ , where $P_d^{(i)}\left(\mathbf{r}\right)$ is the incident power density and $\alpha$ is the fraction of power passing through the skin surface (with the remainder being reflected). In exposure tests, $P_d\left(\mathbf{r}\right)$ is not directly measurable.

• $\mu$ is the skin’s absorption coefficient for the electromagnetic frequency as the wave propagates in the skin. $\mu$ is distinct from the coefficient $\alpha$ introduced above. $\alpha$ represents the fraction of power passing through the skin surface, while $\mu$ refers to the fraction of power absorbed per unit depth as the beam propagates into the skin. $\mu$ has the physical dimension of $1/[\mathrm{length}]$ , whereas $\alpha$ is dimensionless.

• $p_v(\mathbf{r},z)$ is the power absorbed per unit volume at position $(\mathbf{r},z)$ , providing the heating power per unit volume for the temperature evolution. $p_v(\mathbf{r},z)$ decays exponentially in the depth direction and is expressed as:

  $$p_v(\mathbf{r},z)=P_d\left(\mathbf{r}\right)\mu e^{-\mu z}$$

  (1)

• The energy balance gives the governing equation for $T(x,y,z,t)$ :

  $$\underset{\begin{array}{l}\mathrm{rate}\mathrm{of}\mathrm{change}\\ \mathrm{of}\mathrm{heat}\mathrm{in}\mathrm{skin}\end{array}}{\underbrace{{\rho }_m{C}_p\frac{\partial T}{\partial t}}}=\underset{\begin{array}{l}\mathrm{net}\mathrm{heat}\mathrm{influx}\\ \mathrm{from}\mathrm{conduction}\end{array}}{\underbrace{k\left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)}}+\underset{\mathrm{absorbed}\mathrm{power}}{\underbrace{P_d\left(\mathbf{r}\right)\mu e^{-\mu z}}}$$

  (2)

• Effect of lateral heat conduction. In most exposure tests, the length scale in the lateral dimensions is much larger than the penetration depth of the electromagnetic frequency into the skin tissue. In this situation, to the leading-term approximation, lateral heat conduction is negligible. The governing equation becomes:

  $${\rho }_m{C}_p\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial z^2}+P_d\left(\mathbf{r}\right)\mu e^{-\mu z}$$

  (3)

• $T_{\mathrm{base}}$ is the baseline temperature. We assume that the skin temperature is uniform in $(x,y,z)$ at the start of beam exposure ( $t=0$ ), which is referred to as the baseline temperature and is denoted by $T_{\mathrm{base}}$ . This serves as the initial condition for the temperature evolution: $T(x,y,z,0)=T_{\mathrm{base}}$ .

• Zero heat flux at the skin surface. At the skin surface, body heat is slowly lost to the environment through black body radiation and convection cooling by the surrounding air. In homeostasis, this heat loss is balanced by the heat influx carried by blood perfusion. However, for a high-power MMW exposure, the heat loss at the surface is negligible compared to the large energy influx (the power passing into the skin). In our study, we neglect the slow heat loss at the skin surface. For the same reason, the effect of blood perfusion is also neglected in the energy balance of (3). Zero heat flux at the skin surface leads to the boundary condition ${\partial T/\partial z|}_{z=0}=0$ .

• Semi-infinite skin tissue. Given the exponential decay of $p_v(\mathbf{r},z)$ in Equation (1), the power is virtually all absorbed within a few multiples of $1/\mu$ , which is the penetration depth of electromagnetic wave into the skin tissue. In our study, $1/\mu =0.16\mathrm{mm}$ . For the purpose of modeling the temperature evolution, we can approximate the skin tissue as a semi-infinite body, mathematically extending to $z=+\infty$ . The temperature is governed by the initial boundary value problem (IBVP):

  $$\{\begin{array}{c}{\rho }_m{C}_p\frac{\partial T}{\partial t}=k\frac{{\partial }^{2}T}{\partial z^2}+P_d\left(\mathbf{r}\right)\mu e^{-\mu z}\hfill \\ \multicolumn{1}{c}{{\frac{\partial T(x,y,z,t)}{\partial z}|}_{z=0}=0}\\ \multicolumn{1}{c}{T(x,y,z,0)=T_{\mathrm{base}}}\end{array}$$

  (4)

• $T_{\mathrm{act}}$ is the activation temperature of thermal nociceptors.

• $v_{\mathrm{act}}\left(t\right)$ is the activated skin volume at time $t$ . It has the expression:

  $$v_{\mathrm{act}}\left(t\right)\equiv \mathrm{volume}\{(x,y,z)|T(x,y,z,t)\ge T_{\mathrm{act}}\}$$

  (5)

here, $v_{\mathrm{act}}\left(t\right)$ has dependence on $P_d\left(\mathbf{r}\right)$ via $T(x,y,z,t)$ .

• $v_c$ is the threshold on the activated skin volume for flight initiation. In this study, $v_c$ is unknown but is deterministic (i.e., it does not fluctuate from one test to another).

• $t_c$ is the time-of-flight initiation when the activated skin volume reaches the threshold $v_c$ . Given $v_c$ , $t_c$ is solved from the equation:

$${v_{\mathrm{act}}\left(t\right)|}_{t=t_c}=v_c$$

(6)

here, $t_c$ has dependence on $P_d\left(\mathbf{r}\right)$ via $v_{\mathrm{act}}\left(t\right)$ .

• $t_R$ is the human reaction time, also referred to as the SRT.

• $t_F$ is the time-of-flight action, which is the sum of $t_c$ and $t_R$ :

  $$t_F=t_c+t_R$$

  (7)

here, $t_F$ depends on $P_d\left(\mathbf{r}\right)$ through $t_c$ . Although the flight is initiated at time $t_c$ , the occurrence of the flight action is delayed by the human reaction time $t_R$ . The electromagnetic heating continues until $t_F$ , at which point the test subject moves away from the beam.

• The Arrhenius equation describes the reaction rate constant for thermal damage at location $(x,y,z)$ at time $t$ , based on the local temperature $T(x,y,z,t)$ :

$$k_d(x,y,z,t)=A\exp \left(\frac{-\mathrm{\Delta }{E}_a}{RT(x,y,z,t)}\right)$$

The integral of the thermal damage rate constant over time $t$ gives the dimensionless thermal damage parameter, $\mathrm{\Omega }$ , which indicates the degree of burn injury. Since the region near the beam center on the skin surface has the highest temperature, we focus on this region when assessing the skin thermal damage:

$$\mathrm{\Omega }={\int }_0^{t_{\mathrm{stl}}}A\exp \left(\frac{-\mathrm{\Delta }{E}_a}{RT\left(t\right)}\right)\mathrm{dt}$$

(8)

here, $T\left(t\right)$ is the skin surface temperature at the beam center at time $t$ in Kelvin (K). In the integral, the lower limit is the starting time of the beam power, and the upper limit, $t_{\mathrm{stl}}$ , is the settlement time when the skin temperature settles back close to the baseline temperature after the beam power is off. The thermal damage parameter $\mathrm{\Omega }$ depends on $P_d\left(\mathbf{r}\right)$ through $T\left(t\right)$ .

• Parameters in the Arrhenius equation. The parameter $A$ is the frequency factor representing molecular collisions for exciting protein denaturation; $\mathrm{\Delta }{E}_a$ is the activation energy barrier for protein denaturation; and $R$ is the ideal gas constant.

3. Parameter values and normalized formulation

3.1. Parameter values and derived quantities

In the following, we list the parameter values collected from published literature, along with their sources. Some of these parameter values are the same as those used in the ADT CHEETEH model: ²⁰

• Mass density of the skin (dermis): ${\rho }_m=1116\mathrm{kg}/m^3$ .^21,22

• Specific heat capacity of the skin (dermis): $C_p=3300J/(\mathrm{kg}K)$ .^22,23

• Heat conductivity of the skin (dermis): $k=0.445W/(mK)$ .^22,24

• For a Gaussian beam, the power density has the form:

$$P_d\left(\mathbf{r}\right)+P_{\mathrm{d}}^{(0)}\exp \left(\frac{-2|\mathbf{r}{|}^2}{r_b^2}\right)$$

(9)

where $r_b$ is the beam radius, and $P_d^{\left(0\right)}$ is the power density at the beam center absorbed into the skin. $P_d^{\left(0\right)}$ is not directly measurable in exposure tests.

• Penetration depth of a 94-GHz MMW beam into the skin: $1/\mu =0.16\mathrm{mm}$ . ¹³

• Baseline temperature: $T_{\mathrm{base}}=32\mathrm{deg}C=305.15K$ .^20,25

• $T(x,y,z,t)$ , the 3D temperature distribution of the skin at time $t$ , is a derived quantity. In this study, a derived quantity means a quantity calculated using the model. $T(x,y,z,t)$ is solved from IBVP (4) and depends on $P_d\left(\mathbf{r}\right)$ .

• Activation temperature of thermal nociceptors: $T_{\mathrm{act}}=40.4\mathrm{deg}\mathrm{C}=313.55\mathrm{K}$ . ²⁶

• $v_{\mathrm{act}}\left(t\right)$ , the activated skin volume at time $t$ , is a derived quantity. It is calculated from $T(x,y,z,t)$ using formula (5). $v_{\mathrm{act}}\left(t\right)$ has dependence on $P_d\left(\mathbf{r}\right)$ via $T(x,y,z,t)$ .

• $v_c$ , the threshold for the activated skin volume for flight initiation, is unknown.

• $t_c$ , the time-of-flight initiation, is a derived quantity. It is solved in equation $v_{\mathrm{act}}\left(t_c\right)=v_c$ . $t_c$ has dependence on $P_d\left(\mathbf{r}\right)$ via $v_{\mathrm{act}}\left(t\right)$ .

• SRT: $t_R=275\mathrm{ms}$ .^27,28 The mean SRT for women (281 ms) is slightly longer than that for men (268 ms).

• $t_F=t_c+t_R$ , the time-of-flight action, is a derived quantity. It represents the time when the subject starts to move and depends on $P_d\left(\mathbf{r}\right)$ via $t_c$ .

• Parameters in the Arrhenius equation: The frequency factor is $A=8.82\times {10}^{94}{s}^{-1}$ ; the activation energy barrier for protein denature is $\mathrm{\Delta }{E}_a=6.03\times {10}^{5}J/\mathrm{mol}$ ; the ideal gas constant is $R=8.314J/(\mathrm{mol}K)$ .^17,18 It has been argued that the value of the frequency factor $A$ above is unphysically large.^17,18

• $\mathrm{\Omega }$ , the dimensionless thermal damage parameter, is a derived quantity. $\mathrm{\Omega }$ is calculated using formula (8) from the skin temperature $\{T\left(t\right)\}$ over the time course from the start of beam power to the settlement time $t_{\mathrm{stl}}$ , which is when the skin temperature settles back close to the baseline temperature after the beam power is off.

• Classifications of burn injury. Burn injuries are classified into three main categories: first-, second-, and third-degree burns, depending on how deep and how severely they affect the skin tissue. First-degree (or superficial) burns only affect the outer layers of the skin without causing long-term damage. Second-degree (or partial-thickness) burns penetrate deeper than first-degree burns and usually take longer to heal. Third-degree (or full-thickness) burns can cause very serious injuries and require extensive rehabilitation. We adopt the standard classification below, based on the magnitude of the thermal damage parameter $\mathrm{\Omega }$ :^18,23,29

$$\{\begin{array}{c}0.53\le \mathrm{\Omega }<1\to \mathrm{first}-\mathrm{degree}\mathrm{burn}\hfill \\ \multicolumn{1}{c}{1\le \mathrm{\Omega }<{10}^4\to \mathrm{second}-\mathrm{degree}\mathrm{burn}}\\ \multicolumn{1}{c}{\mathrm{\Omega }\ge {10}^4\to \mathrm{third}-\mathrm{degree}\mathrm{burn}}\end{array}$$

(10)

3.2. Normalized formulation

Among all parameters listed in section 3.1, two are unknown or not directly measurable in exposure tests: (1) the volume threshold for flight initiation ( $v_c$ ) and (2) the power density at the beam center absorbed into the skin ( $P_d^{\left(0\right)}$ ). We use the occurrence time-of-flight action $t_F$ to parameterize $P_d^{\left(0\right)}$ . $t_F$ is observed in each exposure test. We normalize the formulation so that the resulting normalized system is independent of the unknown $v_c$ . We select scales for normalization as follows. The penetration depth of a 94-GHz MMW beam into skin provides a length scale in the depth direction, which we denote by $z_s$ :

$$z_s\equiv 1/\mu =0.16\mathrm{mm}$$

To eliminate the unknown volume threshold $v_c$ in the normalized formulation, we select the lateral length scale $r_s$ as:

$$r_s\equiv \sqrt{\frac{v_c}{\pi z_s}}=\sqrt{\frac{\mu v_c}{\pi }}=\mathrm{unknown}$$

Geometrically, $r_s$ is the radius of the cylinder with volume equal to $v_c$ and height equal to $z_s$ . In this way, after normalization, the non-dimensional volume threshold becomes $\pi$ .

The time scale is given by the characteristic time for heat diffusion over a distance of $z_s$ in the depth direction, since lateral heat diffusion is negligible:

$$t_s\equiv \frac{{\rho }_m{C}_p}{k}{z}_s^2=\frac{{\rho }_m{C}_p}{k{\mu }^2}=0.2119s$$

(11)

In the model formulation, we normalize the variables $(\mathbf{r},z,t)$ using the scales $(r_s,z_s,t_s)$ :

$$\{\begin{array}{c}{z}_s\equiv 1/\mu =0.16\mathrm{mm},z_n=\frac{z}{z_s}=\mu z\hfill \\ \multicolumn{1}{c}{t_s\equiv \frac{{\rho }_m{C}_p}{k{\mu }^2}=0.2119s,t_n=\frac{t}{t_s}}\\ \multicolumn{1}{c}{r_s\equiv \sqrt{\frac{\mu v_c}{\pi }},{\mathbf{r}}_n=\frac{\mathbf{r}}{r_s},r_{b,n}=\frac{r_b}{r_s},v_{c,n}=\frac{v_c}{z_s{r}_s^2}=\pi }\end{array}$$

(12)

The temperature is not scaled. We use $T^{(n)}$ to denote the temperature as a function of the scaled variables $(x_n,y_n,z_n,t_n)$ :

$$T^{(n)}(x_n,y_n,z_n,t_n)\equiv T(x,y,z,t)$$

After normalization, the governing equation for $T^{(n)}$ takes the same form as IBVP (4), with the coefficients updated accordingly,

$$\{\begin{array}{c}\frac{\partial T^{(n)}}{\partial t_n}=\frac{{\partial }^2{T}^{(n)}}{\partial z_n^2}+P_d^{(n)}\left({\mathbf{r}}_n\right)e^{-z_n}\hfill \\ \multicolumn{1}{c}{{\frac{\partial T^{(n)}(x_n,y_n,z_n,t_n)}{\partial z_n}|}_{z_n=0}=0}\\ \multicolumn{1}{c}{T^{(n)}(x_n,y_n,z_n,0)=T_{\mathrm{base}}}\end{array}$$

(13)

In Equation (13), the normalized power density has the expression:

$$P_d^{(n)}({\mathbf{r}}_n)=\frac{1}{k\mu }{P}_d(r)=\frac{P_d^{(0)}}{k\mu }\exp \frac{-2|{\mathbf{r}}_{\mathrm{n}}{|}^2}{r_{b,n}^2}$$

(14)

Note that both $T^{(n)}$ and $P_d^{(n)}$ have the physical dimension of temperature, while the normalized beam radius $r_{b,n}$ is dimensionless. Given the parameters $(T_{base},P_d^{\left(0\right)},k,\mu ,r_{b,n})$ , we solve Equation (13) with Equation (14) to calculate the temperature distribution $T^{(n)}(x_n,y_n,z_n,t_n)$ . We then determine the activated volume in the normalized variables $(x_n,y_n,z_n,t_n)$ ,

$$v_{\mathrm{act},n}\left(t_n\right)\equiv \mathrm{volume}\{(x_n,y_n,z_n)|T^{(n)}(x_n,y_n,z_n,t_n)\ge T_{\mathrm{act}}\}$$

(15)

The activated volume in $(x,y,z,t)$ and in $(x_n,y_n,z_n,t_n)$ are related by:

$$v_{\mathrm{act}}\left(t\right)=\left(r_s^2{z}_s\right)v_{\mathrm{act},n}\left(t_n\right)=\frac{v_c}{\pi }{v}_{\mathrm{act},n}\left(t_n\right)$$

(16)

Recall that flight initiation is governed by ${v_{\mathrm{act}}\left(t\right)|}_{t=t_c}=v_c$ , which becomes:

$${v_{\mathrm{act},n}\left(t_n\right)|}_{t_n=t_{c,n}}=\pi$$

(17)

Given the normalized activated volume $v_{\mathrm{act},n}\left(t_n\right)$ , Equation (17) becomes independent of the unknown volume threshold $v_c$ . The effect of the unknown $v_c$ is incorporated into $v_{\mathrm{act},n}\left(t_n\right)$ through the normalized beam radius $r_{b,n}\equiv \frac{r_b}{r_s}$ , where the lateral length scale $r_s$ contains $v_c$ . We specify the beam radius as a multiple of $r_s$ (i.e., by specifying the dimensionless $r_{b,n}$ ), rather than providing a value with physical units. With this approach, the entire normalized formulation is independent of the unknown $v_c$ .

In exposure tests, the beam power is terminated at the time-of-flight action, which is delayed by the human reaction time $t_R$ from the flight initiation. To assess the thermal damage parameter in (8), we need to track the skin temperature from the start of the beam power until the settlement time $t_{\mathrm{stl}}$ , when the skin temperature drops back close to the baseline temperature. In the normalized formulation, we work with the normalized time $t_n$ and $T^{(n)}$ as functions of $t_n$ . The normalized time is related to the physical time by the time scale $t_s$ , such that $t_n=\frac{t}{t_s}$ . We use a change of variables to express $\mathrm{\Omega }$ in terms of the function $T^{(n)}\left(t_n\right)$ :

$$\mathrm{\Omega }=t_s{\int }_0^{t_{\mathrm{stl},n}}A\exp \left(\frac{-\mathrm{\Delta }{E}_a}{R{T}^{(n)}\left(t_n\right)}\right)d{t}_n$$

(18)

The degree of burn injury is determined from $\mathrm{\Omega }$ according to Equation (10).

In summary, we represent the beam radius as a multiple of the unknown length scale $r_s$ , which is derived from the unknown volume threshold $v_c$ . The resulting normalized formulation is independent of the unknown $v_c$ , enabling us to run simulations to assess burn injury risk without requiring knowledge of the value of the volume threshold $v_c$ .

3.3. Parameterization of thermal damage risk

In the previous section, we discussed the advantage of working in the normalized formulation: all parameters in Equation (13) are known and listed in section 3.1, except for the specifications of the Gaussian beam, $(r_b,P_d^{\left(0\right)})$ .

To eliminate the unknown volume threshold $v_c$ in the normalized formulation, the physical beam radius $r_b$ is specified as a multiple of the unknown $r_s$ , which is derived from the unknown $v_c$ . Simulations are carried out with $r_{b,n}=\frac{r_b}{r_s}$ in the normalized formulation. By scaling the lateral length with the unknown $r_s$ , we can compute the thermal damage parameter $\mathrm{\Omega }$ and assess the degree of burn injury without needing the value of the volume threshold $v_c$ .

As described at the end of section 2.1, one approach to parameterizing the thermal damage risk is to use the absorbed power density into the skin, $P_d^{\left(0\right)}$ . Although $P_d^{\left(0\right)}$ is not directly measured in exposure tests, it can be estimated from the specified incident power density and an empirical value for the absorption fraction. The second approach is to use the occurrence time-of-flight action, $t_F$ , to parameterize the thermal damage risk. In real exposure tests, $t_F$ is directly observed as the time when the test subject moves away from the beam. In the absence of randomness in model parameters, there is a one-to-one correspondence between $P_d^{\left(0\right)}$ and $t_F$ , so the two approaches are equivalent in the current deterministic model. However, in the presence of randomness in model parameters, in both approaches, given a fixed value of $P_d^{\left(0\right)}$ (or $t_F$ ), the corresponding thermal damage will be affected by the randomness. In this study, most simulation results are presented in the form of $\mathrm{\Omega }$ versus $t_F$ to emphasize the trade-off between minimizing $t_F$ and avoiding thermal injury. In addition, some results are shown as $\mathrm{\Omega }$ versus $P_d^{\left(0\right)}$ to illustrate the effect of power density on thermal injury.

In Equation (5), the activated volume $v_{\mathrm{act}}\left(t\right)$ increases with time $t$ at any fixed power density $P_d^{\left(0\right)}$ and increases with $P_d^{\left(0\right)}$ at any fixed $t$ . As a result, in Equation (6), a higher value of $P_d^{\left(0\right)}$ implies that the volume threshold $v_c$ is reached at a smaller time-of-flight initiation $t_c$ . Therefore, as a function, $t_c=g(P_d^{\left(0\right)})$ is monotonically decreasing and hence invertible. This monotonic trend allows $t_F$ to be used instead of $P_d^{\left(0\right)}$ to parameterize the thermal damage risk. For a given observed $t_F$ , the underlying $P_d^{\left(0\right)}$ is theoretically determined by the inverse function:

$$P_d^{\left(0\right)}=g^{-1}(t_F-t_R)$$

(19)

In our study, we calculate the thermal damage parameter $\mathrm{\Omega }$ as a function of $(t_F,r_b/r_s)$ . With this setup, the absorbed power density at the beam center, $P_d^{\left(0\right)}$ , is incorporated into $t_F$ , while the effect of the unknown volume threshold $v_c$ is encapsulated in the unknown lateral scale $r_s$ . The unknown $v_c$ is eliminated in the normalized formulation using the length scale $r_s$ , which is derived from $v_c$ . As a result, the physical beam radius is specified as a multiple of the unknown $r_s$ .

3.4. Analytical solution of $T^{(n)}(x_n,y_n,z_n,t_n)$

We begin with the case where the beam power starts at $t_n=0$ and is kept on indefinitely. The temperature as a function of the normalized variables is given by:

$$T^{(n)}(x_{\mathrm{n}},y_{\mathrm{n}},z_{\mathrm{n}},t_{\mathrm{n}})=T_{\mathrm{base}}+\frac{P_d^{(0)}}{k\mu }\mathrm{exp(}\frac{-2|r_n{|}^2}{r_{b,n}^2})U(z_{\mathrm{n}},t_{\mathrm{n}})$$

(20)

where $U(z,t)$ satisfies a parameter-free IBVP of the heat equation:

$$\{\begin{array}{c}\frac{\partial U}{\partial t}=\frac{{\partial }^{2}U}{\partial z^2}+e^{-z}\hfill \\ \multicolumn{1}{c}{{\frac{\partial U(z,t)}{\partial z}|}_{z=0}=0}\\ \multicolumn{1}{c}{U(z,0)=0}\end{array}$$

(21)

IBVP (21) has a closed-form analytical solution:

$$U\left(z,t\right)=\{\begin{array}{cc}-e^{-z}+\frac{e^{-z+t}}{2}erfc\left(\frac{-z+2t}{\sqrt{4t}}\right)+\frac{e^{z+t}}{2}erfc\left(\frac{z+2t}{\sqrt{4t}}\right)-zerfc\left(\frac{z}{\sqrt{4t}}\right)+\frac{2\sqrt{t}}{\sqrt{\pi }}{e}^{\frac{-z^2}{4t}},& t>0\\ 0,& t\le 0\end{array}$$

(22)

In the case where the beam power is terminated at $t_{F,n}$ (the normalized time-of-flight action), the temperature distribution is given by:

$$\begin{array}{c}\multicolumn{1}{c}{T^{(n)}(x_{\mathrm{n}},y_{\mathrm{n}},z_{\mathrm{n}},t_{\mathrm{n}})=T_{\mathrm{base}}+\exp (\frac{-2|r_n^2|}{r_{\mathrm{b},\mathrm{n}}^2}}\end{array})U(z_{\mathrm{n}},t_{\mathrm{n}})-U(z_{\mathrm{n}},t_{\mathrm{n}}-t_{F,n},)$$

(23)

4. Simulations and results

4.1. Fully non-dimensional formulation for simulations

In our simulations, we use the fully non-dimensional formulation. The non-dimensional temperature $T_{\mathrm{nd}}$ and the non-dimensional power density at the beam center $P_{d,\mathrm{nd}}^{\left(0\right)}$ are:

$$\{\begin{array}{c}{T}_{\mathrm{nd}}\equiv \frac{T^{(n)}-T_{\mathrm{base}}}{T_{\mathrm{act}}-T_{\mathrm{base}}},T_{\mathrm{base},\mathrm{nd}}=0\hfill \\ \multicolumn{1}{c}{P_{d,\mathrm{nd}}^{\left(0\right)}\equiv \frac{P_d^{\left(0\right)}}{P_s},P_s\equiv k\mu (T_{\mathrm{act}}-T_{\mathrm{base}})}\end{array}$$

(24)

where $z_n$ , $t_n$ , and ${\mathbf{r}}_n$ are already non-dimensional, as defined in Equation (12). While the beam power is on, the non-dimensional temperature has the expression:

$${\mathrm{T}}_{\mathrm{nd}}(x_{\mathrm{n}},y_{\mathrm{n}},z_{\mathrm{n}},t_{\mathrm{n}})=P_{\mathrm{d},\mathrm{nd}}^{(0)}\exp \left(\frac{-2|{\mathbf{r}}_{\mathrm{n}}{|}^2}{r_{b,n}^2}\right)U(z_{\mathrm{n}},t_{\mathrm{n}})$$

(25)

here, $P_{d,\mathrm{nd}}^{\left(0\right)}$ and $T_{\mathrm{nd}}$ are the non-dimensional, as defined in Equation (24). In terms of $T_{\mathrm{nd}}(x_n,y_n,z_n,t_n)$ , the non-dimensional activated volume, the non-dimensional time-of-flight initiation, and the non-dimensional time-of-flight action are calculated as:

$$\{\begin{array}{c}{v}_{\mathrm{act},n}\left(t_n\right)=\mathrm{vol}\{(x_n,y_n,z_n)|T_{\mathrm{nd}}(x_n,y_n,z_n,t_n)\ge T_{\mathrm{act},\mathrm{nd}}\},T_{\mathrm{act},\mathrm{nd}=1}\hfill \\ \multicolumn{1}{c}{{v_{\mathrm{act},n}\left(t_n\right)|}_{t_n=t_{c,n}}=v_{c,n}{v}_{c,n=\pi }}\\ \multicolumn{1}{c}{t_{F,n}=t_{c,n}+t_{R,n},t_{R,n}=\frac{t_R}{t_s}=\frac{0.275s}{0.2119s}=1.298}\end{array}$$

In exposure tests, the beam power is terminated at the time-of-flight action $t_{F,n}$ , not at the time-of-flight initiation. To correctly assess the thermal damage parameter $\mathrm{\Omega }$ , we track the temperature evolution well beyond the end of the beam power, all the way until the temperature drops back close to the baseline temperature. Over the entire time course, the non-dimensional skin surface temperature at the beam center, $T_{\mathrm{nd}}\left(t_n\right)\equiv T_{\mathrm{nd}}(0,0,0,t_n)$ , follows a unified expression given by:

$$\begin{array}{c}{T}_{\mathrm{nd}}\left(t_n\right)=P_{d,\mathrm{nd}}^{\left(0\right)}(h\left(t_n\right)-h(t_n-t_{F,n})),\mathrm{for}\mathrm{any}{t}_n\hfill \\ \multicolumn{1}{c}{h\left(s\right)\equiv U(0,s)=\{\begin{array}{cc}2\frac{\sqrt{s}}{\sqrt{\pi }}-1+\mathrm{erfcx}\left(\sqrt{s}\right),\hfill & s>0\hfill \\ \multicolumn{1}{c}{0,}& s\le 0\hfill \end{array}}\end{array}$$

where $\mathrm{erfcx}\left(s\right)\equiv e^{s^2}\mathrm{erfc}\left(s\right)$ is the scaled complementary error function. Substituting $T^{\left(n\right)}\left(t_n\right)=T_{\mathrm{base}}+(T_{\mathrm{act}}-T_{\mathrm{base}})T_{\mathrm{nd}}\left(t_n\right)$ into Equation (18), we write the thermal damage parameter $\mathrm{\Omega }$ in terms of the non-dimensional $T_{\mathrm{nd}}\left(t_n\right)$ ,

$$\mathrm{\Omega }=t_s\underset{0}{\overset{t_{\mathrm{stl},n}}{\int }}A\exp \left(\frac{-\mathrm{\Delta }{E}_a}{R[T_{\mathrm{base}}+(T_{\mathrm{act}}-T_{\mathrm{base}})T_{\mathrm{nd}}\left(t_n\right)]}\right)d{t}_n$$

(26)

In numerical integration, Equation (26) is simply an integral of a function of $t_n$ ,

$$\{\begin{array}{c}\mathrm{\Omega }=\underset{0}{\overset{t_{\mathrm{stl},n}}{\int }}{k}_d\left(t_n\right)d{t}_n,k_d\left(t_n\right)\equiv c_1\exp \left(\frac{-1}{c_2+c_3{T}_{\mathrm{nd}}\left(t_n\right)}\right)\hfill \\ \multicolumn{1}{c}{c_1\equiv t_{s}A=\frac{{\rho }_m{C}_p}{k{\mu }^2}A=1.869\times {10}^{94}}\\ \multicolumn{1}{c}{c_2\equiv \frac{R}{\mathrm{\Delta }{E}_a}{T}_{\mathrm{base}}=4.207\times {10}^{-3},}\\ \multicolumn{1}{c}{c_3\equiv \frac{R}{\mathrm{\Delta }{E}_a}(T_{\mathrm{act}}-T_{\mathrm{base}})=1.158\times {10}^{-4}}\end{array}$$

(27)

All quantities in Equation (27) are non-dimensional, including the coefficients $c_1$ , $c_2$ , and $c_3$ . For larger $t_n$ (long after the end of beam power), we observe the following asymptotic behaviors:

$$\begin{array}{c}h\left(t_n\right)\approx \frac{2\sqrt{t_n}}{\sqrt{\pi }}+\cdots \hfill \\ \multicolumn{1}{c}{T_{\mathrm{nd}}\left(t_n\right)\approx \frac{P_{d,\mathrm{nd}}^{\left(0\right)}{t}_{F,n}}{\sqrt{\pi }\sqrt{t_n-t_{F,n}/2}}+\cdots \mathrm{for}{t}_n>t_{F,n}}\end{array}$$

We use the asymptotic formula to estimate the settlement time, $t_{\mathrm{stl},n}$ :

$$\begin{array}{c}{T}_{\mathrm{nd}}\left(t_{\mathrm{stl},n}\right)=c_{\mathrm{stl},n}\equiv 0.5\Rightarrow \frac{P_{d,\mathrm{nd}}^{\left(0\right)}{t}_{F,n}}{\sqrt{\pi }\sqrt{t_n-t_{F,n}/2}}\approx c_{\mathrm{stl},n}\hfill \\ \multicolumn{1}{c}{\Rightarrow t_{\mathrm{stl},n}\approx \frac{1}{2}{t}_{F,n}+{\left(\frac{P_{d,\mathrm{nd}}^{\left(0\right)}{t}_{F,n}}{c_{\mathrm{stl},n}\sqrt{\pi }}\right)}^2}\end{array}$$

(28)

where $c_{\mathrm{stl},n}$ is the coefficient defining the settlement time. Setting $c_{\mathrm{stl},n}=0.5$ gives the time at which the temperature falls back to $(T_{\mathrm{base}}+T_{\mathrm{act}})/2$ .

We use Simpson’s rule to compute the integral in Equation (27). Since the beam power is terminated at $t_{F,n}$ , this discontinuity in the heat source introduces a weak singularity in the temperature at $t_{F,n}$ (the time derivative is discontinuous). To minimize the discretization error in the numerical integration, we set $t_{F,n}$ as a grid point. We use a uniform grid in the interval $[0,t_{F,n}]$ and a non-uniform grid (with gradually increasing grid size) in the interval $[t_{F,n},t_{\mathrm{stl},n}]$ . The union of these two grids covers the entire interval $[0,t_{\mathrm{stl},n}]$ . The integrand function is infinitely differentiable within each interval. Specifically, we define the two grids as:

$$\{\begin{array}{c}{t}_j^{\left(1\right)}=\frac{j}{N_1}{t}_{F,n},j=0,1,2,\dots ,N_1\hfill \\ \multicolumn{1}{c}{t_j^{\left(2\right)}=t_{F,n}{(1-(1-\beta )\frac{j}{N_2})}^{-2},j=0,1,2,\dots ,N_2,}\\ \multicolumn{1}{c}{\beta \equiv \sqrt{\frac{t_{F,n}}{t_{\mathrm{stl},n}}}}\end{array}$$

(29)

The thermal damage parameter $\mathrm{\Omega }$ , as given in Equation (27), is discretized as:

$$\begin{array}{l}\mathrm{\Omega }=\underset{0}{\overset{t_{F,n}}{\int }}{k}_d\left(t_n\right)d{t}_n+\underset{t_{F,n}}{\overset{t_{stl,n}}{\int }}{k}_d\left(t_n\right)d{t}_n\approx {\mathrm{\Omega }}_1+{\mathrm{\Omega }}_2\\ {\mathrm{\Omega }}_1=\sum _{j=0}^{N_1-1}(k_d\left(t_j^{\left(1\right)}\right)+4k_d\left(t_{j+\frac{1}{2}}^{\left(1\right)}\right)+k_d\left(t_{j+1}^{\left(1\right)}\right))\frac{{\left(\Delta t\right)}_j^{\left(1\right)}}{6}\\ t_{j+\frac{1}{2}}^{\left(1\right)}\equiv \frac{1}{2}\left(t_j^{\left(1\right)}+t_{j+1}^{\left(1\right)}\right),{\left(\Delta t\right)}_j^{\left(1\right)}\equiv t_{j+1}^{\left(1\right)}-t_j^{\left(1\right)}\\ {\mathrm{\Omega }}_2=\sum _{j=0}^{N_2-1}(k_d\left(t_j^{\left(2\right)}\right)+4k_d\left(t_{j+\frac{1}{2}}^{\left(2\right)}\right)+k_d\left(t_{j+1}^{\left(2\right)}\right))\frac{{\left(\Delta t\right)}_j^{\left(2\right)}}{6}\\ t_{j+\frac{1}{2}}^{\left(2\right)}\equiv \frac{1}{2}\left(t_j^{\left(2\right)}+t_{j+1}^{\left(2\right)}\right),{\left(\Delta t\right)}_j^{\left(2\right)}\equiv t_{j+1}^{\left(2\right)}-t_j^{\left(2\right)}\end{array}$$

4.2. Simulation results based on the listed parameters

Given the parameter values listed in section 3.1, an exposure test in our model is completely specified by $(t_F,r_b/r_s)$ . In this setting, the absorbed power density at the beam center, $P_d^{\left(0\right)}$ , is encapsulated in $t_F$ , while the effect of the unknown volume threshold, $v_c$ , is embedded in the unknown lateral scale, $r_s$ . This setup is essential for bypassing the unknowns $(P_d^{\left(0\right)},v_c)$ .

We first demonstrate how the beam radius $r_b/r_s$ influences $T\left(t\right)$ , the skin surface physical temperature at the beam center as a function of physical time $t$ .

In Figure 1, the observed time-of-flight action is fixed at $t_F=1\mathrm{second}$ . The left panel displays the results of an exposure test with $r_b=r_s$ , where $r_s$ represents the unknown lateral length scale derived from the unknown volume threshold $v_c$ . The skin surface temperature increases monotonically over the interval $[0,t_F]$ , driven by electromagnetic heating. After the beam power is terminated at $t_F$ , the temperature decreases steadily for $t>t_F$ . The maximum temperature reached is $T_{\max }=63.39\mathrm{deg}\mathrm{C}$ . The temperature remains relatively high for a substantial period of time after $t_F$ . The corresponding absorbed power density is $P_d^{\left(0\right)}=5.17W/c{m}^2$ , and the thermal damage parameter is $\mathrm{\Omega }=2.505$ , indicating a second-degree burn ( $\mathrm{\Omega }>1$ ). The right panel of Figure 1 presents the results for a slightly larger beam radius, $r_b=1.25r_s$ . The maximum temperature of the time course decreases to $T_{\max }=56.06\mathrm{deg}\mathrm{C}$ , significantly lower than that shown in the left panel for $r_b=r_s$ . The absorbed power density drops to $P_d^{\left(0\right)}=3.96W/c{m}^2$ , and the thermal damage parameter decreases to $\mathrm{\Omega }=0.026$ , which is well below the threshold for a first-degree burn ( $\mathrm{\Omega }>0.53$ ). Comparing the two panels of Figure 1, we observe that the same flight action time of $t_F=1s$ can be achieved with a significantly lower thermal injury risk by using a larger beam radius with a lower power density. A larger beam radius heats a larger skin area and reduces the activated depth required to reach the volume threshold. Here, the activated depth refers to the depth of the activated 3D skin region. Since the electromagnetic heating source decays exponentially with depth, a smaller activated depth necessitates a lower power density, which results in a lower maximum temperature and, consequently, a reduced thermal damage parameter $\mathrm{\Omega }$ .

Figure 1.

Results for $t_F=1s$ . Skin surface temperature versus $t$ . (a) $r_b=r_s$ ; (b) $r_b=1.25r_s$ .

In Figure 2, for $r_b=r_s$ , we examine how the thermal damage parameter $\mathrm{\Omega }$ and the absorbed power density $P_d^{\left(0\right)}$ vary with the observed time-of-flight action $t_F$ . In real exposure tests, the time-of-flight action $t_F$ is not precisely prescribed or controlled in the experimental setup. Instead, $t_F$ is observed, reflecting the underlying absorbed power density $P_d^{\left(0\right)}$ . The experimenter can adjust the power density at the beam source to influence $P_d^{\left(0\right)}$ , and, consequently, to influence $t_F$ . As the observed $t_F$ increases, both $\mathrm{\Omega }$ and $P_d^{\left(0\right)}$ decrease monotonically. At a fixed beam radius, the activated depth required to reach the volume threshold is also fixed. However, for a larger observed $t_F$ (and thus a larger $t_c$ ), heat conduction has more time to smooth out the temperature distribution along the depth. A more uniform temperature along the depth means that a given activated depth is achieved with a lower maximum temperature, which in turn reduces the thermal damage parameter $\mathrm{\Omega }$ .

Figure 2.

Results for $r_b=r_s$ . (a) Thermal damage parameter $\mathrm{\Omega }$ versus observed flight action time $t_F$ ; (b) Absorbed power density $P_d^{\left(0\right)}$ versus $t_F$ .

Next, we investigate the effect of beam radius $r_b$ on $\mathrm{\Omega }$ versus $t_F$ and on $\mathrm{\Omega }$ versus $P_d^{\left(0\right)}$ . Figure 3 shows results for four values of $r_b$ , ranging from $r_b=0.5r_s$ to $r_b=20r_s$ .

Figure 3.

$\mathrm{\Omega }$ versus $t_F$ at $r_b=0.5r_s$ , $r_s$ , $2r_s$ , and $20r_s$ .

Similar to what was demonstrated in Figure 2, at each fixed $r_b$ , the thermal damage parameter $\mathrm{\Omega }$ is always a decreasing function of $t_F$ . In all panels, the vertical axis for $\mathrm{\Omega }$ uses the same logarithmic scale with thresholds labeled for first-, second-, and third-degree burns. The horizontal axis for $t_F$ is set differently for each individual $r_b$ to capture the transition of $\mathrm{\Omega }$ from no injury (below the first-degree burn) to second-degree burn and above. As the beam radius $r_b$ increases, the transition location in $t_F$ decreases rapidly. For $r_b=0.5r_s$ , the transition spans $t_F=33~40s$ (the top left panel of Figure 3). With a beam of radius of $r_b=0.5r_s$ , it is impossible for flight action to occur below 30 s while avoiding thermal injury. When the beam radius is increased to $r_b=r_s$ , the transition occurs within a narrow interval around $t_F\approx 1.15s$ (the top right panel of Figure 3). When the beam radius is further increased to $r_b=2r_s$ , the transition location drops further to $t_F\approx 0.45s$ (the bottom left panel of Figure 3). It is clear that if we want flight action to occur at a small $t_F$ while avoiding thermal injury, a large beam radius should be used. There is, however, a lower bound on $t_F$ . Since $t_F=t_c+t_R$ and $t_R=0.275s$ is independent of the beam specifications, the smallest achievable $t_F$ is definitely bounded from below by $0.275s$ . This obvious lower bound, however, is not attainable when we add the requirement of avoiding thermal injury. To make $t_F\approx t_R$ , we need to push the flight initiation time $t_c=t_F-t_R$ to zero, which requires a very large absorbed power density to heat the skin from its baseline temperature ( $T_{\mathrm{base}}$ ) to the nociceptor activation temperature ( $T_{\mathrm{act}}$ ) almost instantly. This very large absorbed power density does not end at the flight initiation $t_c$ . It continues until the flight action at $t_F=t_c+t_R$ . This very large absorbed power density over a time period of $t_F>0.275s$ will cause significant thermal injury. For the same reason, when $t_F$ is not far above $t_R$ , a small decrease in $t_F$ results in a large percentage decrease in $t_c=t_F-t_R$ , which requires a large percentage increase in the absorbed power density $P_d^{\left(0\right)}$ , leading to a significant increase in the thermal damage parameter $\mathrm{\Omega }$ . That is why the curve of $\mathrm{\Omega }$ against $t_F$ is nearly vertical for $t_F$ approaching $t_R$ , as shown in the bottom right panel of Figure 3. With the requirement of avoiding thermal injury, $t_F$ approaching $t_R$ is not achievable, regardless of the beam radius and power density used. A practically achievable lower bound on $t_F$ that avoids thermal injury is about $t_F\approx 0.347s$ , which requires a fairly large beam radius of $r_b=20r_s$ .

In Figure 4, we plot $\mathrm{\Omega }$ versus $P_d^{\left(0\right)}$ for the same set of $r_b$ values used in Figure 3. For each $r_b$ , the thermal damage parameter $\mathrm{\Omega }$ increases monotonically with $P_d^{\left(0\right)}$ , the realized heating power density (the absorbed power density). As $P_d^{\left(0\right)}$ increases, there is a transition from no injury to injury. The transition location in $P_d^{\left(0\right)}$ increases monotonically with beam radius $r_b$ . The increase is more pronounced at small $r_b$ . The transition location increases from $P_d^{\left(0\right)}=0.5W/c{m}^2$ at $r_b=0.5$ to $P_d^{\left(0\right)}=4.2W/c{m}^2$ at $r_b=1$ . Note that in Figure 4, $\mathrm{\Omega }$ varies gradually with $P_d^{\left(0\right)}$ , unlike the steep curve of $\mathrm{\Omega }$ against $t_F$ for $t_F$ approaching $t_R$ , observed in Figure 3. For example, for $r_b=20r_s$ in the bottom right panel of Figure 4, a first-degree burn occurs at the absorbed power density $P_d^{\left(0\right)}\approx 10.5(W/c{m}^2)$ , while a third-degree burn won’t happen until $P_d^{\left(0\right)}\approx 17.5(W/c{m}^2)$ .

Figure 4.

$\mathrm{\Omega }$ versus $P_d^{\left(0\right)}$ at $r_b=1.5r_s$ , $2r_s$ , $5r_s$ , and $20r_s$ .

In Figure 5, we compare six curves of $\mathrm{\Omega }$ versus $t_F$ for beam radii ranging from $r_b=0.75r_s$ to $r_b=20r_s$ , all plotted in a single panel.

Figure 5.

$\mathrm{\Omega }$ versus $t_F$ for several values of beam radius $r_b$ .

As the beam size $r_b$ is increased from slightly below $r_s$ to $r_s$ , and then to slightly above $r_s$ , the location of the transition window from injury to no injury drops very rapidly in $t_F$ . When $r_b=0.75r_s$ , $t_F=3s$ is definitely not safe for avoiding thermal injury. When $r_b=1.25r_s$ , it is safe even at $t_F=1s$ . As $r_b$ increases further above $r_s$ , the location of the transition window converges to a practical lower bound of about $t_F=0.347s$ . This attainable lower bound is notably close to the human reaction time $t_R=0.275s$ . As discussed above, the attainable lower bound on $t_F$ is not only limited by $t_R$ , but also constrained by the requirement to avoid thermal injury.

4.3. Sensitivity study

We study the sensitivity of $\mathrm{\Omega }$ versus $t_F$ when the parameters are perturbed from their listed values in section 3.1. In our sensitivity study, we vary one parameter at a time while keeping the others fixed at their listed values.

Recall that in our model formulation, the beam is completely described by $(t_F,r_b/r_s)$ , where the absorbed power density $P_d^{\left(0\right)}$ is implicit in the observed time-of-flight action $t_F$ . The beam radius is specified as a multiple of the lateral length scale $r_s\equiv \sqrt{\frac{\mu v_c}{\pi }}$ , which varies with the parameters $\mu$ and $v_c$ , but is independent of other parameters. As a result, the sensitivity study with respect to $\mu$ or $v_c$ requires additional adjustments when setting the dimensionless parameters. In the sensitivity study below, the physical beam radius is fixed when a parameter is perturbed. We first examine sensitivity with respect to the parameters that do not affect the lateral length scale $r_s$ . We set $r_b=r_s$ .

We start with the human reaction time $t_R$ . Curves of $\mathrm{\Omega }$ versus $t_F$ for values of $t_R$ ranging from 0.138 s to 0.55 s are plotted in Figure 6. The listed value of $t_R$ in section 3.1 is $t_R=0.275s$ . As shown in Figure 6, at any given time-of-flight action $t_F$ , the thermal damage parameter $\mathrm{\Omega }$ increases monotonically with the human reaction time $t_R$ . This result is intuitive. At a fixed $t_F=t_c+t_R$ , a larger $t_R$ gives a smaller time to flight initiation $t_c$ , which requires a higher power density for more rapid heating. A higher heating power over a fixed heating duration $t_F$ increases the thermal damage risk.

Figure 6.

$\mathrm{\Omega }$ versus $t_F$ for several values of human reaction time $t_R$ . $r_b=r_s$ .

We explore the effect of the skin baseline temperature $T_{\mathrm{base}}$ . The listed value of $T_{\mathrm{base}}$ in section 3.1 is $T_{\mathrm{base}}=32\mathrm{deg}C$ . Curves of $\mathrm{\Omega }$ versus $t_F$ for values of $T_{\mathrm{base}}$ ranging from $29\mathrm{deg}\mathrm{C}$ to $35\mathrm{deg}\mathrm{C}$ are plotted in Figure 7, which shows that at any given time-of-flight action $t_F$ , the thermal damage parameter $\mathrm{\Omega }$ decreases as the skin baseline temperature $T_{\mathrm{base}}$ increases. A higher skin baseline temperature $T_{\mathrm{base}}$ implies that a smaller temperature increase is needed to reach the nociceptor activation temperature. At a fixed $t_c=t_F-t_R$ , a smaller temperature increase translates to a lower power density needed, which in turn decreases the thermal damage risk, since the heating duration $t_F$ is fixed.

Figure 7.

$\mathrm{\Omega }$ versus $t_F$ for several values of $T_{\mathrm{base}}$ . $r_b=r_s$ .

We continue by examining the effect of the nociceptor activation temperature $T_{\mathrm{act}}$ . The listed value of $T_{\mathrm{act}}$ in section 3.1 is $T_{\mathrm{act}}=40.4\mathrm{deg}\mathrm{C}$ . Figure 8 shows curves of $\mathrm{\Omega }$ versus $t_F$ for values of $T_{\mathrm{act}}$ ranging from $37.4\mathrm{deg}\mathrm{C}$ to $43.4\mathrm{deg}\mathrm{C}$ . At any given time-of-flight action $t_F$ , the thermal damage parameter $\mathrm{\Omega }$ increases monotonically with the nociceptor activation temperature $T_{\mathrm{act}}$ . Similar to the situation with $T_{\mathrm{base}}$ , a higher $T_{\mathrm{act}}$ means a larger temperature increase is needed to reach $T_{\mathrm{act}}$ from $T_{\mathrm{base}}$ , which translates to a higher power density and thus a higher thermal damage risk.

Figure 8.

$\mathrm{\Omega }$ versus $t_F$ for several values of nociceptor activation temperature $T_{\mathrm{act}}$ . $r_b=r_s$ .

Next, we investigate the effect of the skin heat conductivity $k$ . The listed value of $k$ in section 3.1 is $k=0.445W/(\mathrm{mK})$ . We plot curves of $\mathrm{\Omega }$ versus $t_F$ for values of $k$ ranging from $0.223W/(\mathrm{mK})$ to $0.89W/(\mathrm{mK})$ in Figure 9. At any given time-of-flight action $t_F$ , the thermal damage parameter $\mathrm{\Omega }$ decreases as the skin heat conductivity $k$ increases. We need to interpret this result within our model framework. At a fixed $t_c=t_F-t_R$ , a larger conductivity $k$ more effectively smooths the temperature distribution in depth via heat conduction, making the temperature increase along the depth more uniform. A more uniform heating lowers the power density needed to reach the required activated depth and reduces the maximum skin temperature, which ultimately decreases the thermal damage risk.

Figure 9.

$\mathrm{\Omega }$ versus $t_F$ for several values of skin heat conductivity $k$ . $r_b=r_s$ .

Another quantity that influences the effectiveness of heat conduction in smoothing out the temperature in depth is the volumetric heat capacity $\rho C_p$ . In our model, the parameters $\rho$ and $C_p$ appear only in the combination $\rho C_p$ . The listed value of $\rho C_p$ in section 3.1 is $\rho C_p=3.68\times {10}^{6}J/(m^{3}K)=3.68J/({\mathrm{cm}}^{3}K)$ . Figure 10 shows the curves of $\mathrm{\Omega }$ versus $t_F$ for values of $\rho C_p$ ranging from $1.84J/({\mathrm{cm}}^{3}K)$ to $7.37J/({\mathrm{cm}}^{3}K)$ . At any given time-of-flight action $t_F$ , the thermal damage parameter $\mathrm{\Omega }$ increases monotonically with the volumetric heat capacity $\rho C_p$ . The mechanism behind this behavior is similar to that for conductivity $k$ in Figure 9. At a fixed conductivity $k$ , a larger volumetric heat capacity $\rho C_p$ renders the heat conduction less effective in smoothing out the temperature in depth because more heat needs to be transferred to reduce a given temperature gradient. Thus, a larger volumetric heat capacity $\rho C_p$ has an effect similar to that of a smaller conductivity $k$ .

Figure 10.

$\mathrm{\Omega }$ versus $t_F$ for several values of volumetric heat capacity $\left(\rho C_p\right)$ . $r_b=r_s$ .

Finally, we examine the effect of the electromagnetic penetration depth $1/\mu$ and the unknown volume threshold $v_c$ . In the model formulation, the beam radius is specified as a multiple of the lateral length scale $r_s\equiv \sqrt{\frac{\mu v_c}{\pi }}$ , which depends on the parameters $\mu$ and $v_c$ . In the sensitivity study, to isolate the effect of $\mu$ , we need to fix all other parameters, including the physical beam radius. Operationally, while we vary $\mu$ in simulations, it is desirable to continue working with the normalized formulation described in section 3.2. To accommodate both of these, we use the listed value $1/{\mu }^{\left(0\right)}=0.16\mathrm{mm}$ to define a reference lateral length scale $r_s^{\left(0\right)}\equiv \sqrt{\frac{{\mu }^{\left(0\right)}{v}_c}{\pi }}$ , which remains unaffected when $\mu$ is varied in the sensitivity study. The reference scale $r_s^{\left(0\right)}$ is solely for the purpose of specifying $r_b$ . In the sensitivity study, we fix the physical beam radius at $r_b=r_s^{\left(0\right)}$ . In simulations based on the normalized formulation, as $\mu$ is varied, the normalized beam radius $r_{b,n}$ is adjusted to reflect the underlying fixed physical beam radius:

$$r_{b,n}\equiv \frac{r_b}{r_s}=\frac{r_s^{\left(0\right)}}{r_s}\frac{r_b}{r_s^{\left(0\right)}}=\sqrt{\frac{{\mu }^{\left(0\right)}}{\mu }}\mathrm{when}{r}_b\mathrm{is}\mathrm{fixed}\mathrm{at}{r}_s^{\left(0\right)}$$

(30)

In Equation (30), the physical beam radius is fixed at $r_b=r_s^{\left(0\right)}$ . When $1/\mu$ is increased by a factor of 2, the dimensionless $r_{b,n}$ increases by a factor of $\sqrt{2}$ . In the normalized formulation, $\mu$ also affects the time scale $t_s=\frac{{\rho }_m{C}_p}{k{\mu }^2}$ , the power density scale $P_s\equiv k\mu (T_{\mathrm{act}}-T_{\mathrm{base}})$ and the coefficient $c_1=\frac{{\rho }_m{C}_p}{k{\mu }^2}A$ in (27) for evaluating $\mathrm{\Omega }$ .

Curves of $\mathrm{\Omega }$ versus $t_F$ for the electromagnetic penetration depth $1/\mu$ ranging from 0.08 mm to 0.32 mm are plotted in Figure 11. The listed value in section 3.1 is $1/\mu =0.16\mathrm{mm}$ . As shown in Figure 11, at any given time-of-flight action $t_F$ , the thermal damage parameter $\mathrm{\Omega }$ decreases as $1/\mu$ is increased. A larger $1/\mu$ implies a more uniform heat source in the depth. When heating the skin to a given activated depth, a more uniform heating leads to a more uniform temperature and a lower maximum skin temperature, which decreases the thermal damage parameter $\mathrm{\Omega }$ .

Figure 11.

$\mathrm{\Omega }$ versus $t_F$ for several values of penetration depth $1/\mu$ . $r_b=r_s^{\left(0\right)}$ .

In the sensitivity study with respect to the unknown volume threshold $v_c$ , we fix the physical beam radius as follows. We select a reference $v_c^{\left(0\right)}$ to define a reference lateral length scale $r_s^{\left(0\right)}\equiv \sqrt{\frac{\mu v_c^{\left(0\right)}}{\pi }}$ . Since $v_c$ is unknown, the reference $v_c^{\left(0\right)}$ may be arbitrary. We investigate the effect of varying $v_c$ relative to this reference value. Again, the reference length scale $r_s^{\left(0\right)}$ is solely for the purpose of specifying the physical beam radius $r_b$ . In the sensitivity study, we fix the physical beam radius at $r_b=r_s^{\left(0\right)}$ . In simulations based on the normalized formulation, as $v_c$ is varied, the normalized beam radius $r_{b,n}$ is adjusted to reflect the same underlying fixed physical beam radius:

$$r_{b,n}\equiv \frac{r_b}{r_s}=\frac{r_s^{\left(0\right)}}{r_s}\frac{r_b}{r_s^{\left(0\right)}}=\sqrt{\frac{v_c^{\left(0\right)}}{v_c}}\mathrm{when}{r}_b\mathrm{is}\mathrm{fixed}\mathrm{at}{r}_s^{\left(0\right)}$$

(31)

In Equation (31), the physical beam radius is fixed at $r_b=r_s^{\left(0\right)}$ . When $v_c$ is increased by a factor of 2, the dimensionless $r_{b,n}$ decreases by a factor of $\sqrt{2}$ . The most prominent difference between the effects of $\mu$ and $v_c$ is that $v_c$ affects only the lateral length scale $r_s$ , not other parameters. It follows that the effect of increasing $v_c$ while fixing the physical beam size $r_b$ is the same as that of reducing $r_b$ while fixing $v_c$ . Figure 12 shows the curves of $\mathrm{\Omega }$ versus $t_F$ for the volume threshold $v_c$ ranging from $0.5v_c^{\left(0\right)}$ to $2v_c^{\left(0\right)}$ , where $v_c^{\left(0\right)}$ is the reference value. At any given time-of-flight action $t_F$ , the thermal damage parameter $\mathrm{\Omega }$ increases monotonically with $v_c$ . The trend of $\mathrm{\Omega }$ increasing with $v_c$ in Figure 12 corresponds to the trend of $\mathrm{\Omega }$ decreasing with $r_s$ observed earlier in Figure 5.

Figure 12.

$\mathrm{\Omega }$ versus $t_F$ for several values of volume threshold $v_c$ . $r_b=r_s^{\left(0\right)}$ .

5. Conclusion

In this study, we formulated a model to assess the skin thermal damage risk during exposure tests, where the electromagnetic beam radiation persists until the occurrence of flight action, at which point the test subject moves away from the beam. The increase in skin temperature is driven by the absorbed electromagnetic power and is governed by an IBVP of the heat equation. Whenever the local temperature exceeds the activation threshold, thermal nociceptors in the skin are activated, sending a nerve signal. The burning sensation corresponds to the aggregated nerve signals from the exposed skin spot. Once the activated skin volume reaches a certain threshold, the nociceptor signal becomes strong enough to trigger flight initiation. However, there is a time delay before the flight action occurs, as it takes time for the nociceptor signal to travel to the brain, for the flight signal to be processed in the brain, for it to propagate to the muscles, and for the muscles to actuate the flight. The delay between flight initiation and observed flight action is known as the human reaction time. During this period, the beam power continues to be applied until the flight action occurs.

In the thermal model, a relevant beam specification is the power density absorbed by the skin, $P_d^{\left(0\right)}$ , which is not directly measurable in exposure tests. However, the time-of-flight action, $t_F$ , is reliably observed during these tests. In our model, $t_F$ is a strictly decreasing function of $P_d^{\left(0\right)}$ . This monotonic relationship allows us to use the observable $t_F$ as a surrogate for $P_d^{\left(0\right)}$ to parameterize the thermal injury risk.

We collected parameter values from the published literature for the thermally induced flight model. The only exception is the threshold $v_c$ , which represents the activated skin volume required for flight initiation, and remains unknown. To eliminate the unknown volume threshold $v_c$ in the normalized formulation, we normalized the independent variables. In the depth direction, the length scale $z_s$ is naturally determined by the penetration depth of the electromagnetic frequency. The time scale $t_s$ is derived from heat diffusion and the depth scale, as heat diffusion is significant only in the depth direction. The lateral length scale $r_s$ is derived from the unknown volume threshold $v_c$ such that after the scaling, the normalized volume threshold is set to $\pi$ (the volume of a cylinder of height 1 and radius 1). Since $v_c$ is unknown, $r_s$ is also unknown. However, when the physical beam radius is specified as a multiple of the unknown lateral scale $r_s$ , there are no unknown parameters in the normalized formulation. This setup enables simulations in the normalized formulation to compute the thermal damage parameter $\mathrm{\Omega }$ and use it to classify the burns: first-degree burn ( $\mathrm{\Omega }\ge 0.52$ ), second-degree burn ( $\mathrm{\Omega }\ge 1$ ), and third-degree burn ( $\mathrm{\Omega }\ge {10}^4$ ).

In simulations of thermal injury risk, $\mathrm{\Omega }$ is parameterized by the observed time-of-flight action $t_F$ . The heating power density $P_d^{\left(0\right)}$ from the beam is embedded in the observed $t_F$ . This aligns well with the realistic scenario in exposure tests, where $t_F$ is measured, but $P_d^{\left(0\right)}$ is not directly observed. Our simulations lead to the following conclusions:

• At a fixed physical beam radius, a longer observed time-of-flight action corresponds to a smaller underlying heating power density. This results in more uniform heating in depth, as the effect of heat conduction becomes more pronounced. Consequently, when the activated volume threshold is reached, the maximum temperature is lower, and the thermal damage parameter is reduced.

• At a fixed observed time-of-flight action, using a larger beam radius in an exposure test heats a larger skin area, reducing the required activated depth to reach flight initiation. As a result, the underlying heating power density needed to reach flight initiation is smaller, which leads to a lower thermal damage parameter, since the heating duration remains fixed.

• At a fixed observed time-of-flight action, when the human reaction time is shorter, the time to flight initiation becomes longer. This results in a lower underlying heating power density needed to reach flight initiation, which in turn leads to a smaller thermal damage parameter, since the total heating duration remains fixed.

• At a fixed observed time-of-flight action, when the skin baseline temperature is higher or the nociceptor activation temperature is lower, the temperature difference between the baseline and the activation threshold becomes smaller. As a result, the underlying heating power density required to reach flight initiation decreases, leading to a smaller thermal damage parameter, since the heating duration is fixed.

• At a fixed observed time-of-flight action, when the skin heat conductivity is higher or the volumetric heat capacity is lower, heat diffusion becomes more effective in smoothing out the temperature gradient in depth. As a result, the activated depth required for flight initiation is reached with a lower maximum skin temperature, which yields a smaller thermal damage parameter.

• At a fixed observed time-of-flight action, when the penetration depth of the electromagnetic frequency into the skin is larger (i.e., when the skin absorption coefficient for the electromagnetic frequency is smaller), the heating in the depth becomes more uniform. This results in the activated depth needed for flight initiation being reached with a lower maximum skin temperature, which reduces the thermal damage parameter.

• At a fixed observed time-of-flight action, when the threshold for the activated skin volume required for flight initiation is smaller, the heating power density needed to reach this volume threshold is lower. As a result, the thermal damage parameter becomes smaller since the heating duration remains fixed.

• Although a smaller volume threshold and a larger penetration depth both reduce the thermal damage parameter, they operate through different mechanisms and have varying magnitudes of effect. A smaller volume threshold directly reduces the heating power density required for flight initiation, making it a more immediate and effective way to lower the thermal damage parameter. In contrast, a larger penetration depth increases the uniformity of electromagnetic heating in the depth direction. While it does not directly reduce the heating power density required, it lowers the thermal damage parameter by reducing the maximum skin temperature, resulting in a smaller temperature gradient in the depth. However, the effect of a larger penetration depth on lowering the thermal damage parameter is less pronounced.

Finally, it is important to highlight the limitations of the model in the current study. The current mathematical formulation relies on the following assumptions:

1. The skin tissue is modeled as a single layer with uniform material properties, particularly assuming they are independent of depth within the skin.

2. All thermal nociceptors in the skin are assumed to behave identically, including (1) responding solely to the instantaneous temperature, (2) having the same deterministic activation temperature, and (3) transducing an electrical signal of the same magnitude when activated by a local temperature that exceeds the activation threshold.

3. Thermal nociceptors are uniformly distributed throughout the skin, such that the sum of nociceptor signals is directly proportional to the activated skin volume.

4. The heat sensation is proportional to the sum of nociceptor signals, which in turn is proportional to the activated skin volume. When the activated skin volume exceeds a subjective threshold, a flight signal is triggered.

5. There is a time gap between the initiation of the flight signal and the actual flight action (i.e., when the subject begins to move). This time gap is akin to the SRT measured in psychophysical experiments.

6. Once the subject starts to move, the beam radiation is considered terminated.

7. The thermal damage is modeled using the Arrhenius equation.

Even when assumption (1) holds, the parameter values found in the literature carry significant uncertainty. The effect of this uncertainty and fluctuations in model parameters is investigated in our sensitivity study, where we quantify how the injury risk varies when a parameter is perturbed. Assumptions (2) to (7) are included primarily to simplify the current model, enabling it to have an analytical solution with a small set of parameters. This simplicity offers the advantage of operational applicability, even in the absence of detailed knowledge about skin structure and properties.

In addition to the assumptions listed above, our current study uses the Arrhenius equation for thermal damage assessment based on the skin surface temperature. However, it has been shown that the surface temperature correlates poorly with thermal damage from MMW radiation in real exposure tests, especially for third-degree burns. ³⁰ For a more accurate assessment, it is necessary to calculate the thermal damage profile in the depth direction and assign burn injury based on this profile, rather than relying solely on the surface value. ³¹

In our subsequent studies, we will explore incorporating several factors into the model, including (a) multi-layer skin tissue, (b) non-uniform distribution of thermal nociceptors, (c) fluctuations in the thermal nociceptor activation temperature, (d) additional types of thermal nociceptors that respond to heat rate rather than instantaneous temperature, and (e) thermal damage assessment based on the temperature profile in the depth.