**2. Problem formulation**

A two-dimensional squared-shaped biological tissue having the dimensions of 2 mm � 2 mm was considered in the present study, which is shown in **Figure 1**. The short-pulse laser is vertically incident at the top surface of the domain. Short laser pulses have been modeled as Gaussian profile in the space domain and Top-Hat profile in the time domain with a pulse width of *tp*. When light propagates through the biological tissue, the intensity has two components: collimated and diffuse. The collimated component of the intensity has been mathematically expressed by Eq. (6) [21].

$$I\_{\varepsilon}(\mathbf{x}, \mathbf{y}, \boldsymbol{\Omega}, t) = I\_{c\max} \times \exp\left\{- (\mathbf{x} - \mathbf{L}/2)^2 / d^2 \right\} \times \exp\left\{- \beta (\boldsymbol{W} - \boldsymbol{y}) / |\boldsymbol{\eta}\_c|\right\} \tag{6}$$

$$\begin{split} \times & \left[ H\{\beta \varepsilon t - \beta (\boldsymbol{W} - \boldsymbol{y}) / |\boldsymbol{\eta}\_c|\} - H\{ (\beta \varepsilon t - \beta (\boldsymbol{W} - \boldsymbol{y}) / |\boldsymbol{\eta}\_c|) - \beta \varepsilon t\_p \} \right] \\ & \times \delta(\boldsymbol{\mu} - \boldsymbol{\mu}\_c) \delta(\boldsymbol{\eta} - \boldsymbol{\eta}\_c) \end{split}$$

where, *Icmax*, *L*, *W*, *μ*, and *η* represent the maximum intensity incident on the top surface of the tissue phantom, the length of tissue, the width of the tissue, and direction cosines in *x*- and *y*-directions, respectively. Here, *δ* denotes the Dirac-delta function, which considers the angle of incidence, and *H* is the Heaviside function which ensures that the laser energy is available at any location only for the pulse width of the laser. The first and second exponential terms in Eq. (6) represent the Gaussian profile in the space domain and the attenuation of collimated intensity due to absorption and out-scattering, respectively.

### **Figure 1.** *Schematic diagram of the physical domain under consideration (dimensions in mm).*

The diffuse component of the intensity is obtained by solving the RTE. So, Eq. (1) is simplified for the two-dimensional rectangular domain and also neglects the emission from the tissue phantom because the intensity of the short-pulse laser used for irradiating the sample is much higher than the blackbody radiation intensity:

$$\frac{1}{c}\frac{\partial I\_d}{\partial t} + \mu \frac{\partial I\_d}{\partial \mathbf{x}} + \eta \frac{\partial I\_d}{\partial \mathbf{y}} = -\beta I\_d + \frac{\sigma\_t}{4\pi} \int\_{4\pi} \mathbf{I} \otimes (\boldsymbol{\Omega}, \boldsymbol{\Omega}') d\boldsymbol{\Omega}' \tag{7}$$

The boundary conditions for a diffusely emitting and reflecting wall at a given point *rw* on that surface can be expressed as [21]:

$$I(r\_w, \hat{s}, t) = \in (r\_w)I\_b(r\_w, t) + \frac{\rho(r\_w)}{\pi} \int\_{\vec{n} \cdot \hat{s}' < 0} I(r\_w, \hat{s}', t) |\hat{n} \bullet \hat{s}'| d\mathcal{Q}', \hat{n} \bullet \hat{s} > 0 \tag{8}$$

Here *n*^ is the inward normal vector to the surface. The first and second terms on the right-hand side of Eq. (8) represent the emission and reflection from the wall, respectively. The expression for boundary conditions given by Eq. (8), in general, holds valid for all the four walls of the rectangular enclosure.

Once the intensity distribution inside the laser-irradiated biological tissue is obtained, then the divergence of the radiative heat flux (Eq. 5) is calculated. As earlier mentioned, finding the temperature distribution inside the biological tissue subjected to the short-pulse laser is a multi-time scale problem. So, we need to divide the Eq. (2) into Eqs. (9) and (10). Eq. (9) is solved for determining local temperature rise due to a single pulse for which the time scale is the order of picoseconds.

$$
\rho c\_v \frac{\partial T}{\partial t} = -\nabla \bullet \overrightarrow{q\_r} \tag{9}
$$

To determine the temperature distribution inside the biological tissue subjected to a train of pulses, the temperature rise due to the single pulse is added to the temperature distribution at the previous time instant if the time difference between the two consecutive pulses is equal to the time step employed else it is not added. The temperature distribution at the previous time instant is used as the initial condition to determine the temperature distribution at the current time instant by solving the Eq. (10) [21]. It is to be noted that the time scale used for solving Eq. (10) is the order of milliseconds.

$$
\rho c\_v \frac{\partial T}{\partial t} = -\nabla \bullet \overrightarrow{q} + a\_b \rho\_b c\_b (T\_b - T) + Q\_m \tag{10}
$$

The DPL model-based bio-heat transfer equation has been solved in the present study, which is obtained by eliminating temperature from Eq. (10) using the Eq. (4), and the governing equation is written in terms of heat flux which is expressed as:

$$\frac{\tau\_q}{a}\frac{\partial^2 \overrightarrow{q}}{\partial t^2} + \frac{1}{a} \frac{\partial \overrightarrow{q}}{\partial t} = \nabla \left(\nabla \bullet \overrightarrow{q}\right) + \tau\_T \frac{\partial \left[\nabla \left(\nabla \bullet \overrightarrow{q}\right)\right]}{\partial t} \tag{11}$$

It is to be noted here that the blood perfusion and metabolic heat generation terms have been neglected while deriving Eq. (11) because their contribution is negligible to temperature rise [21].

The required boundary conditions for solving the Eq. (10) are as follows: the all three walls except the top wall (shown in **Figure 1**) are maintained at the core body temperature (37°C), while the top wall is subjected to the convective boundary conditions (heat transfer coefficient of 15 W/m2 ∙K and ambient temperature of 25°C) [21]. The initial condition (*t =* 0) for temperature is equal to 37°C, and the two initial conditions are required for Eq. (11) because it has a second-order derivative concerning the time, so heat flux and time derivative of heat flux are assumed to be zero.

Eqs. (10) and (11) have been numerically solved using the FVM, and its solution is given in Section 3.

The following laser parameters have been used in this study: beam diameter (*d*)= 0.025 mm, amplitude of pulse (*Icmax*)= 1.6�10�<sup>3</sup> J/mm<sup>2</sup> /ps, pulse width (*tp*)= 4.6667 ps, repetition rate (*fr*)= 1 kHz and wavelength= 1100 nm [21]. The absorption coefficient (*κ*) and scattering coefficient ð Þ *<sup>σ</sup><sup>s</sup>* of tissue are 0.051 mm�<sup>1</sup> and 6.14 mm�<sup>1</sup> , respectively [4, 21]. The thermo-physical properties of the tissue are as follows: density (*ρ*)= 1000 kg/m<sup>3</sup> , thermal conductivity (*k*)= 0.63 W/m∙K and specific heat (*cv*)= 4200 J/kg∙K [21].
