Thermal Effect of a Revolving Gaussian Beam on Activating Heat-Sensitive Nociceptors in Skin

Abstract

We consider the problem of inducing withdrawal reflex on a test subject by exposing the subject’s skin to an electromagnetic beam. Heat-sensitive nociceptors in the skin are activated wherever the temperature is above the activation temperature. Withdrawal reflex occurs when the activated volume reaches a threshold. Previously we studied static beams with 3 types of power density distribution: Gaussian, super-Gaussian, and flat-top. We found that the flaptop is the best and the Gaussian is the worst in their performance with regard to 1) minimizing the time to withdrawal reflex, 2) minimizing the energy consumption and 3) minimizing the maximum temperature increase. The less-than-desirable performance of Gaussian beams is attributed to the uneven distribution of power density resulting in low energy efficiency: near the beam center the high power density does not contribute proportionally to increasing the activated volume; outside the beam effective radius the low power density fails to activate nociceptors. To overcome the drawbacks of Gaussian beams, in this study, we revolve a Gaussian beam around a fixed point to make the power density more uniformly distributed. We optimize the performance over two parameters: the spot size of static beam and the radius of beam revolution. We find that in comparison with a static Gaussian beam, a revolving Gaussian beam can reduce the energy consumption, and at the same time lower the maximum temperature.

Share and Cite:

Wang, H. , Burgei, W. 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. doi: 10.4236/jamp.2021.91007.

1. Introduction

Common applications of radiofrequency (RF) radiation such as laser surgery and cancer hyperthermia lead to personnel exposure to electromagnetic energy [1]. Due to the shallow penetration depth of RF radiation in biological tissues, the energy and, therefore, heat associated with the electromagnetic wave is deposited in tissues near the surface of skin. We consider the thermal effect of RF radiation on human skin. Thermal responses of humans to RF radiation have been an active research topic for several decades [2].

The electromagnetic energy deposited by RF radiation increases the skin temperature. Heat-sensitive nociceptors in the skin are activated wherever the temperature is above the activation temperature. When the activated volume reaches a threshold withdrawal reflex occurs. We examine 3 aspects of an electromagnetic beam: the time to withdrawal reflex, the energy consumption and the maximum temperature at withdrawal reflex. In this study, we extend our previous analysis of static beams [3] to a Gaussian beam revolving around a fixed point. We find that in comparison with a static Gaussian beam, a revolving Gaussian beam has lower energy consumption and lower maximum temperature at withdrawal reflex.

2. Mathematical Formulation for a Moving Beam

We consider the situation where a skin area of the test subject is exposed to a moving beam [4] [5] [6]. We adopt a formulation similar to the one in our previous studies [3] [7]. Let T ( z , r , t ) denote the temperature of skin as a function of 3-D spatial coordinates ( z , r ) and time t. Here z is the depth from the skin surface and r is the 2-D coordinates on the skin surface. We assume 1) the moving beam is perpendicular to the skin surface (i.e., incident angle = 0); 2) before the exposure, the 3-D skin has a uniform initial temperature T base (the baseline temperature); and iii) heat conduction is included only in the depth direction.

At each surface location r , the temperature distribution T ( z , r , t ) is governed by

{ ρ m C p T t = k 2 T z 2 + P ( r , t ) μ exp ( μ z ) T z | z = 0 = 0, T ( z ,0 ) = T base (1)

where

· ρ 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,

· μ is the absorption coefficient of the skin for the beam frequency, and

· P ( r , t ) is the beam power density absorbed into the skin at location r at time t; for a moving beam whose center changes with time, P ( r , t ) is a function of t.

The electromagnetic energy absorbed into the skin increases the temperature. We adopt the same general assumption used in the ADT CHEETEH model for linking the thermal effect to withdrawal reflex: heat-sensitive nociceptors are activated wherever the local temperature is above the nociceptor activation temperature T act ; when a sufficient number of nociceptors are activated, withdrawal reflex occurs and the test subject moves away from the beam [8] [9]. We consider the case where heat-sensitive nociceptors are uniformly distributed in skin. In this case, the number of nociceptors activated is directly proportional to the volume of the activated region ( T T act ). Withdrawal reflex occurs when the activated volume is above a critical threshold v c . In general, threshold v c has some randomness, varying with the test subject and/or the exposure location on the body. In addition, it may vary with the internal and environmental factors that affect human physiologically or psychologically. This randomness models the underlying uncertainty in the occurrence of withdrawal reflex [10]. Here for simplicity, we consider the case of v c = deterministic .

We first non-dimensionalize variables and functions in (1). The depth scale is provided by, 1 / μ , the characteristic scale of electromagnetic energy penetrating in the depth direction; the length scale for surface coordinates is derived from volume threshold v c .

Length scale in the depth direction and time scale:

z s 1 μ , t s ρ m C p k μ 2

Non-dimensional depth and time:

z nd z z s = μ z , t nd t t s = t k μ 2 ρ m C p

Length scale for surface coordinates and volume scale:

r s v c π z s , v s r s 2 z s = v c π

Non-dimensional surface coordinates and volume:

r nd r r s , V nd V v s

Temperature scale and power density scale:

Δ T T act T base , P s k μ Δ T

Non-dimensional temperature as a function of ( z nd , r nd , t nd ) :

T nd ( z nd , r nd , t nd ) T ( z , r , t ) T base Δ T

Non-dimensional power density as a function of ( r nd , t nd ) :

P nd ( r nd , t nd ) P ( r , t ) P s

In the non-dimensional formulation, model parameters have simple values

v c,nd = π , T base,nd = 0, T act,nd = 1

At each surface location r nd , the non-dimensional temperature is governed by

{ T nd t nd = 2 T nd z nd 2 + P nd ( r nd , t nd ) exp ( z nd ) T nd z nd | z nd = 0 = 0 , T ( z nd ,0 ) = 0 (2)

The solution of initial boundary value problem (2) has the analytical expression

{ T ( z , r , t ) = 0 t P ( r , s ) G ( z , t s ) d s G ( z , t ) = 1 2 erfc ( 2 t z 4 t ) e t z + 1 2 erfc ( 2 t + z 4 t ) e t + z (3)

Here, for conciseness, we have dropped the subscript “nd” and used the simple notations for all non-dimensional quantities. For example, T ( z , r , t ) in (3) means T nd ( z nd , r nd , t nd ) . Based on the non-dimensional temperature, we calculate the non-dimensional activated volume

V act ( t ) = Volume { ( z , r ) | T ( z , r , t ) T act } , T act = 1

By definition, the (non-dimensional) reflex time t ref is the time it takes for the activated volume to reach threshold v c . Using v c = π , we write the equation for t ref as

Volume { ( z , r ) | T ( z , r , t ref ) 1 } = π (4)

3. Temperature Solution for a Gaussian Beam Revolving around a Fixed Point

For a Gaussian beam, the power density distribution is axisymmetric with respect to the beam center. At location r on the skin, the time dependent power density from the moving Gaussian beam has the expression

P ( r , t ) = p GB ( | r r ctr ( t ) | ) (5)

where

· r ctr ( t ) is the (non-dimensional) center of the moving Gaussian beam,

· p GB ( r ) is the power distribution function of the static Gaussian beam,

p GB ( r ) = 2 P tot π r eff 2 exp ( 2 ( r r eff ) 2 ) = P tot r eff 2 p GB ( 0 ) ( r r eff ) (6)

· r eff is the (non-dimensional) effective mode radius of the beam [11], and

· P tot p GB ( r ) 2 π r d r is the (non-dimensional) total power of the beam,

· p GB ( 0 ) ( r ) 2 π exp ( 2 r 2 ) is the power distribution function of the standardized Gaussian beam with effective radius r eff ( 0 ) = 1 and total power P tot ( 0 ) = 1 .

Note that a moving Gaussian beam is completely specified by ( { r ctr ( t ) } , r eff , P tot ) . Non-dimensional r eff and P tot are related to their corresponding physical quantities by

r eff = r eff , phy r s , P tot = P tot , phy P s r s 2

where P tot denotes the total power and P s the scale for power density. We consider a Gaussian beam moving along a circle with uniform speed, as shown in Figure 1. For convenience, we define the center of the circle as the origin, and the initial position of the beam center as the x-axis. Let

· R be the (non-dimensional) radius of the circular trajectory, and

· ω be the (non-dimensional) angular velocity of the beam center.

Consider a location on the skin surface with coordinates r = r ( cos θ , sin θ ) . At time t, the beam center is r ctr ( t ) = ( R cos ( ω t ) , R sin ( ω t ) ) . As described in (5), the power density at location r is determined by the distance to the beam center,

P ( r , t ) = p GB ( d ( t ) )

d ( t ) | r r ctr ( t ) | = r 2 + R 2 2 r R cos ( ω t θ ) (7)

The distance d ( t ) is periodic and even in variable ϕ ( ω t θ ) . As a result, power density p GB ( d ( t ) ) has a Fourier cosine expansion in ϕ

p GB ( d ( t ) ) = a 0 2 + k = 1 a k cos ( k ϕ ) , ϕ ( ω t θ ) (8)

With the Fourier expansion, we write temperature solution (3) as

T ( z , r , t ) = a 0 2 0 t G ( z , t s ) d s + k = 1 a k 0 t cos ( k ( ω s θ ) ) G ( z , t s ) d s I k (9)

Figure 1. Schematic diagram of a Gaussian beam moving along a circle.

Applying integration by parts to I k , we get

I k = 1 k ω 0 t G ( z , t s ) d sin ( k ( ω s θ ) ) = 1 k ω ( sin ( k ( ω t θ ) ) G ( z ,0 ) + sin ( k θ ) G ( z , t ) + 0 t sin ( k ( ω s θ ) ) G ( z , t s ) t d s )

It follows that | I k | O ( 1 / ( k ω ) ) . Fourier coefficients in (8) have the expression

a k = 1 π 0 2 π cos ( k ϕ ) p GB ( r 2 + R 2 2 r R cos ϕ ) d ϕ (10)

Since p GB ( r 2 + R 2 2 r R cos ϕ ) is infinitely differentiable, coefficient a k decays very rapidly with k, satisfying a k = o ( 1 / k M ) for any M 1 . Combining a k = o ( 1 / k M ) and | I k | O ( 1 / ( k ω ) ) , we obtain | k = 1 a k I k | = O ( 1 / ω ) and we write solution (9) as

T ( z , r , t ) = p rev ( | r | ) 0 t G ( z , s ) d s + O ( 1 ω ) (11)

where p rev ( | r | ) is the power density at location r averaged over one revolution of the revolving Gaussian beam. p rev ( r ) has the expression

p rev ( r ) = 1 2 π 0 2 π p GB ( r 2 + R 2 2 r R cos ϕ ) d ϕ = P tot r eff 2 1 2 π 0 2 π p GB ( 0 ) ( ( r r eff ) 2 + ( R r eff ) 2 2 ( r r eff ) ( R r eff ) cos ϕ ) d ϕ (12)

Notice that in (12), p rev ( r ) is proportional to a function of ( r / r eff ), with ( R / r eff ) as a parameter, and is independent of θ . In Figure 2 and Figure 3, we compare the power density p GB ( r ) of a static Gaussian beam and the average power density p rev ( r ) of the corresponding revolving Gaussian beam.

Figure 2. (a) Power distribution of the Gaussian beam at its initial position and (b) Average power distribution of the revolving Gaussian beam.

Figure 3. Comparison of p GB ( r ) and p rev ( r ) . The total beam power is P tot = 4 , the beam radius is r eff = 0.425 and the radius of revolution is R = 0.5 .

It is evident in Figure 3 that moving a Gaussian beam along a properly selected circle is an effective way of spreading out the beam energy to heat a large region more uniformly and at the same time to avoid overheating the center. Note that in Figure 3, the horizontal axis is the radial coordinate; power densities p GB ( r ) and p rev ( r ) are conserved in the sense 2 π r p GB ( r ) d r = 2 π r p rev ( r ) d r .

4. Performance and Optimization of Fast Revolving Gaussian Beams

We consider the case of a fast revolving Gaussian beam (i.e., ω 1 ). Neglecting the O ( 1 / ω ) term in (11), the temperature solution becomes

T ( z , r , t ) = p rev ( | r | ) H ( z , t ) , H ( z , t ) 0 t G ( z , s ) d s (13)

With solution (13), we rewrite the equation for t ref given in (4) as

Volume { ( z , r ) | p rev ( | r | ) H ( z , t ref ) 1 } = π (14)

Here, although not explicitly indicated in the notation, p rev ( | r | ) is affected by the radius of beam revolution (R) and the radius of static beam spot ( r eff ). Below we optimize the performance of a fast revolving Gaussian beam over parameters ( R , r eff ) .

We fix the total power P tot of the beam instead of fixing the peak power density. For a static Gaussian beam, distribution function (6) dictates that the peak power density is inversely proportional to the squares of beam radius: p GB ( 0 ) P tot / r eff 2 . For a revolving Gaussian beam, the total power is further spread out, resulting in a more uniform power density (see Figure 3). We first look at the time to withdrawal reflex vs beam radius for a static Gaussian ( R = 0 ), shown in the left panel of Figure 4. The corresponding maximum temperature increase is plotted in the right panel of Figure 4. When the beam radius is small,

Figure 4. Results of a static beam with a fixed total beam power. (a) Reflex time vs beam radius and (b) Maximum temperature vs beam radius.

the power is too much concentrated in a small area. Since the penetration depth of electromagnetic energy into the skin is limited, heat propagation in the depth direction is carried by conduction, which takes time. As a result, a small beam radius needs a long time to reach the activated volume threshold and leads to large surface temperature increase. When the beam radius is large, the power is spread out too thin over a large area. With a low power density, it takes long time to reach the activation temperature. At any fixed value of P tot , there is an optimal beam radius ( r opt , marked by dash-dotted lines in Figure 4) for minimizing the reflex time. For a static Gaussian beam, minimizing the reflex time produces a fairly large temperature increase: max T nd > 5.0 (see the right panel of Figure 4). Here we are studying the system in the non-dimensional formulation and thus all quantities are dimensionless.

To demonstrate the effect of revolving the beam, we select 3 values of r eff around the optimal beam radius ( r opt ) and examining the reflex time vs the radius of beam revolution in the left panel of Figure 5. The corresponding maximum temperature increase is plotted in the right panel of Figure 5. Results of the static beam with r opt are shown as dash horizontal lines in Figure 5. For r eff > r opt , the beam power is already spread out too thin in the static beam. Any further spreading by revolving the beam leads to a larger reflex time (solid red line in Figure 5). For r eff < r opt , moving the beam along a small circle increases the heating area and reduces the reflex time. When the radius of revolution is large, the beam power is spread out too thin and it takes long time to reach the activation temperature. It follows that for r eff < r opt , there is an optimal radius of revolution ( R opt , marked by dotted lines in Figure 5) for minimizing the reflex time. The minimum reflex time achieved by revolving a beam with r eff < r opt may be lower than what can be achieved with a static beam with r opt (comparing solid blue and dashed lines in Figure 5). Notice that by revolving a beam with r eff < r opt , both a smaller reflex time and a smaller maximum temperature can be achieved simultaneously. For a beam of r eff = 0.5 , minimizing the reflex time by revolving the beam yields a temperature increase of max T nd 2.35 , much lower than that for the static beam with the optimal radius (dashed horizontal line in Figure 5).

Next, we minimize the reflex time over the 2D parameter space ( r eff , R ) at each given value of P tot . The left panel of Figure 6 shows the optimal ( r eff , R ) vs P tot . The ratio of optimal r eff to optimal R is plotted in the right panel. The optimal r eff and R increase with the total power. For a larger total beam power, the reflex time is smaller and the heat conduction in the depth direction has less time to take its effect. As a result, the heating in the depth direction is less effective and the power needs to be spread out to heat a larger area in order to reach the activated volume threshold in the shortest time. We point out that the ratio of optimal r eff to optimal R stays roughly unchanged (≈0.85) over a wide range

Figure 5. Results of a revolving beam with a fixed total beam power. (a) Reflex time vs radius of revolution and (b) Maximum temperature vs radius of revolution.

Figure 6. Optimization of a revolving beam over the 2D parameter space ( r eff , R ) . (a) Optimal ( r eff , R ) vs P tot and (b) (Optimal r eff )/(optimal R) vs P tot .

of P tot . When we need to increase the spreading of beam power, to make the resulting power density as uniform as possible, it is best to increase beam radius r eff and radius of beam revolution R proportionally. Recall that in the expression of power distribution for a revolving beam given in (12), p rev ( r ) is proportional to a function of ( r / r eff ) with ( R / r eff ) as a parameter. When the ratio ( r eff / R ) is fixed, power distribution p rev ( r ) for various combinations of ( r eff , R ) is congruent to each other via horizontal/vertical linear scaling. Figure 6 tells us that the ratio ( r eff / R ) 0.85 gives the optimal shape of power distribution. The blue solid line in Figure 3 shows this optimal shape.

In Figure 7, we compare the performance of revolving Gaussian beam with those of static Gaussian, static super-Gaussian, and static flat-top beams. All beams are optimized with respect to r eff (and R, for the revolving beam only). In the left panel, the total energy consumption of revolving Gaussian is lower than that of static Gaussian but higher than those of static super-Gaussian and static flat-top. In the right panel, the maximum temperature of revolving Gaussian is lower than that of static Gaussian but higher than those of static super-Gaussian and static flat-top. Figure 7 demonstrates that by moving a Gaussian beam along a circle, we can lower both the energy consumption and the maximum temperature. The improvement is fairly uniform over a wide range of

total beam power, from P tot = 1 2 to P tot = 64 (non-dimensional).

In Figure 6 and Figure 7, the revolving Gaussian beam has been optimized over parameter space ( r eff , R ) to achieve the minimum reflex time (and equivalently, the minimum energy consumption). Given the Gaussian form of beam power distribution, it is not possible to further reduce the reflex time/energy consumption. However, it is possible to reduce the maximum temperature at the price of a slightly larger energy consumption. As demonstrated in Figure 4 and

Figure 7. Comparison of a revolving Gaussian beam with 3 static beams: Gaussian, super-Gaussian and flat-top. (a) P tot vs total energy consumption at optimal ( r eff , R ) and (b) P tot vs maximum temperature at optimal ( r eff , R ) .

Figure 8. Performance of the revolving Gaussian beam with ( r eff , R ) = β × optimal ( r eff , R ) . (a) Total energy consumption vs P tot and (b) Maximum temperature vs P tot .

Figure 5, the maximum temperature is reduced when the beam power is spread out more by increasing the beam radius ( r eff ) and/or increasing the radius of beam revolution (R). We increase r eff and R proportionally. In Figure 8, we compare the performance of revolving Gaussian beams with ( r eff , R ) = β × optimal ( r eff , R ) for several values of β 1 . As β is increased from 1, the maximum temperature is reduced at the price of increased energy consumption. At β = 1.3 , the revolving Gaussian beam has a maximum temperature of max T nd 2.0 , similar to that of the optimal static flat-top beam which has the lowest maximum temperature among static beams. In particular, the maximum temperature caused by the revolving Gaussian beam with β = 1.3 is much smaller than max T nd > 5.0 of the optimal static Gaussian beam. At the same time, the energy consumption of the revolving Gaussian beam with β = 1.3 is still noticeably lower than that of the optimal static Gaussian beam (comparing dash-dotted curve in Figure 8 and dashed blue curve in Figure 7).

5. Concluding Remarks

We studied the problem of causing a heat-induced withdrawal reflex on a test subject by exposing the subject’s skin to a Gaussian electromagnetic beam. Heat-sensitive nociceptors in the skin are activated wherever the temperature is above the activation temperature. Withdrawal reflex occurs when the activated volume reaches a threshold. The Gaussian beam has a highly uneven distribution of power density, resulting in low energy efficiency: near the beam center the high power density does not contribute proportionally to increasing the activated volume; outside the beam effective radius the power density is not high enough to activate nociceptors and is wasted. In addition, the high power density near beam center produces a large local temperature increase that does not contribute to increasing the activated volume but poses a serious burn injury risk for the test subject. To overcome the drawbacks of static Gaussian beams, we explored the idea of making the power density more uniform by moving a Gaussian beam along a circle. Given the total beam power, the two adjustable parameters of a revolving Gaussian beam are the radius of static beam and the radius of beam revolution. We optimized over these two parameters to minimize the time to withdrawal reflex (i.e. to minimize the energy consumption and maximize the energy efficiency). We found that in comparison with a static Gaussian beam, a revolving Gaussian beam can lower the reflex time and the energy consumption and at the same time reduce the maximum temperature at withdrawal reflex.

Acknowledgements and Disclaimer

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] Ryan, K.L., D’Andrea, J.A. Jauchem, J.R. and Mason, P.A. (2000) Radio Frequency Radiation of Millimeter Wave Length: Potential Occupational Safety Issues Relating to Surface Heating. Health Physics, 78, 170-181.
https://doi.org/10.1097/00004032-200002000-00006
[2] Foster, K.R., Ziskin, M.C. and Balzano, Q. (2016) Thermal Response of Human Skin to Microwave Energy: A Critical Review. Health Physics, 111, 528-541.
https://doi.org/10.1097/HP.0000000000000571
[3] 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
[4] 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
[5] 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
[6] 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
[7] 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
[8] Cazares, S.M., Snyder, J.A., Belanich, J., Biddle, J.C., Buytendyk, A.M., Teng, S.H.M. and O’Connor, K. (2019) Active Denial Technology Computational Human Effects End-to-End Hypermodel for Effectiveness (ADT CHEETEH-E). Institute for Defense Analyses. IDA Document NS D-10435.
https://doi.org/10.1007/s41314-019-0023-7
[9] Teng, S.H.M., Biddle, J.C., Snyder, J.A., Belanich, J., Buytendyk, A.M., O’Connor, K. and Cazares, S.M. (2019) Fidelity Improvements and Additional Analyses for the Active Denial Technology Computational Human Effects End-to-End Hypermodel (ADT CHEETEH). Institute for Defense Analyses. IDA Document D-10762.
[10] Wang, H., Burgei, W.A. and Zhou, H. (2018) Dose-Injury Relation as a Model for Uncertainty Propagation from Input Dose to Target Dose. American Journal of Operations Research, 8, 360-385.
https://doi.org/10.4236/ajor.2018.85021
[11] Paschotta, R. (2008) Article on “Effective Mode Area” in the Encyclopedia of Laser Physics and Technology. Wiley-VCH.

Copyright © 2024 by authors and Scientific Research Publishing Inc.

Creative Commons License

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.