# Numerical Study of Photoacoustic Pressure for Cancer Therapy
## 2 Photoacoustic Model  
This section presents the physical model and the numerical method used to solve the photoacoustic problem. The photoacoustic or thermoacoustic process refers to the generation of acoustic waves by absorption of electromagnetic energy. Such an absorption of the radiation by a nanoparticle acts like a heat source that depends on the space \(\mathbf{x}\) and time \(t\) and is given by:  
\[Q(\mathbf{x},t)=\frac{\omega}{2}\epsilon_{0}Im\left(\epsilon_{r}(\mathbf{x}) \right)|\mathbf{E}(\mathbf{x},t)|^{2},\quad\forall\mathbf{x}\in\Omega, \tag{1}\]  
where \(\omega\) is the angular frequency of the incoming wave, \(\mathbf{E}\) the electric field, \(\epsilon_{0}\) is the permittivity of vacuum, \(Im(\epsilon_{r})\) is the imaginary part of \(\epsilon_{r}\) the relative complex permittivity of materials that are functions of the 3D spatial coordinates \(\mathbf{x}\) in the domaine \(\Omega\). The laser illumination pulse satisfies the electromagnetic equation:  
\[\
abla\times\left(\frac{1}{\epsilon_{r}(\mathbf{x})}\
abla\times\mathbf{E}( \mathbf{x},t)\right)-\frac{1}{c^{2}}\frac{\partial^{2}\mathbf{E}(\mathbf{x},t )}{\partial t^{2}}=\mathbf{0},\quad\forall\mathbf{x}\in\Omega, \tag{2}\]  
where (\(\
abla\times\)) is the curl operator and \(c\) the velocity of light. Equation (2) is solved with a radiation boundary condition [11; 12] at the external border \(\Gamma\) of the computationnal domain \(\Omega\) corresponding to the free propagation of the scattered field and the incident laser pulse of amplitude:  
\[E_{0}=\left[\frac{2c\mu_{0}E_{laser}}{\pi R_{beam}^{2}\tau_{pulse}}\right]^{1/ 2}, \tag{3}\]  
where \(\mu_{0}\) is the permeability of the vacuum, \(E_{laser}\) is the laser energy, \(R_{beam}\) is the radius of the laser beam and \(\tau_{pulse}\) the duration of the laser pulse. From such a heat source, the temperature elevation is governed by the equation:  
\[\left[\rho(\mathbf{x})C_{p}(\mathbf{x})\frac{\partial}{\partial t}-\
abla. \left(\kappa(\mathbf{x})\
abla\right)\right]T(\mathbf{x},t)=Q(\mathbf{x},t), \quad\forall\mathbf{x}\in\Omega, \tag{4}\]  
where \((\
abla. )\) and \(\
abla\) are the divergence and the gradient operators, \(\rho\) is the mass density of the material, \(C_{p}\) is the specific heat capacity, \(\kappa\) the thermal conductivity and \(T\) is the temperature. These are functions of time and of the 3D spatial coordinates \(\mathbf{x}\). Therefore, the resolution of the thermal problem requires solving the heat equation (Equation (4)) with a source \(Q\) produced by the incoming electromagnetic illumination (Equation (1)) with the boundary condition \(T(\mathbf{x},t)=T_{0}\) a the boundary of the computational domain (i.e., \(\forall\mathbf{x}\in\Gamma\)). Such a thermal expansion induces acoustic displacement and the equation for photoacoustic generation is derived [2]:  
\[\left[\frac{1}{\rho(\mathbf{x})v^{2}(\mathbf{x})}\frac{\partial^{2}}{\partial t ^{2}}-\
abla.\left(\frac{1}{\rho(\mathbf{x})}\
abla\right)\right]P(\mathbf{x},t)=\beta(\mathbf{x})\frac{\partial^{2}}{\partial t^{2}}T(\mathbf{x},t),\quad \forall\mathbf{x}\in\Omega, \tag{5}\]  
with \(v\) the sound velocity in the medium, \(\beta\) is the expansion coefficient and \(P\) is the acoustic pressure with initial and boundary conditions:  
\[P(\mathbf{x},t=0)=p_{0}\quad\text{ and }\quad n.\left(\frac{1}{\rho(\mathbf{x})} \
abla\right)P(\mathbf{x},t)=0,\quad\forall\mathbf{x}\in\Gamma, \tag{6}\]  
with \(p_{0}\) being the initial pressure and \(n\) the outgoing normal vector. The resolution of the weak coupled photothermal problem gives the spatio-temporal distribution of the pressure in the computational domain. Equations (1)-(5) are numerically solved with a 3D Finite Element Method (FEM) with adaptive remeshing scheme using an a priori error estimator based on the Hessian of the solution calculated atthe previous step of remeshing.