Absorption of the laser beam by the medium subsequently gives rise to an absorbed optical energy, $\gamma I$. Assuming that all the optical energy is converted to heat, a local increase in temperature occurs. Temperature distribution, $T$, in the absence of convection and phase transition, can be computed using the following heat equation:
\begin{equation}
    k \nabla ^2 T = \rho C \frac{\partial T}{\partial t} - \gamma I
    \label{eq:eqChaleur}
\end{equation}
where $k$ is the thermal conductivity, $\rho$ is the density, $C$ is the heat capacity and $t$ is the time. Calculating the exact solution to this equation is beyond the scope of this article, but we can roughly approximate the first and second terms to be $k T / \gamma^2$ and $\rho C T / \tau$, respectively, during laser emission. Given that $k$ = 0.6 W.m$^{-1}$.K$^{-1}$ (water thermal conductivity), $\rho$ = 1000 kg.m$^{-3}$ (water density), $C$ = 4180 J.kg$^{-1}$.m$^{-3}$ (water heat capacity), $\gamma^{-1} \approx$ 40 $\mu$m and $\tau$ = 10 ns, the first term is negligible compared to the second one. Substituting low-energy experimental parameters ($E$ = 10 mJ, $S$ = 20 mm$^{2}$) leads to a maximum increase in temperature of 3 K, which produces a local dilatation of the medium. The induced displacements can then generate shear waves, which constitutes the \textit{thermoelastic regime}.

To estimate the initial displacement amplitude in this regime, we assumed the medium as homogeneous and isotropic. As the depth of absorption (about 40 $\mu$m) is 100 times smaller than the beam diameter (5 mm), we discarded any boundary effects. The stress, $\sigma_{zz}$, is the sum of the axial strain component and the thermal expansion component \cite{scruby1990laser}:
\begin{equation}
    \sigma_{zz} = (\lambda + 2 \mu) \frac{\partial u_z}{\partial z} - 3(\lambda + \frac{2}{3}\mu) \alpha \frac{ E}{\rho C S \zeta}
    \label{eq:stressThermo}
\end{equation}
where $\lambda$ and $\mu$ are respectively the first and second Lamé's coefficients, $\alpha$ is the thermal dilatation coefficient, and $\zeta$ is the average depth of absorption. This equation can be simplified by the fact that in most soft media, including biological tissues, $\mu \ll \lambda$. Moreover, in the absence of external constraints normal to the surface, the stress across the surface must be zero, i.e. $\sigma_{zz} (z=0) = 0$. This allows the equation \ref{eq:stressThermo} to be integrated, giving at the surface a displacement $u_z = 3 \alpha E/(\rho C S)$. Substituting the same experimental parameters used previously along with $\alpha$ = 70.10$^{-6}$ K$^{-1}$ (water linear thermal dilatation coefficient), we obtain a displacement $u_z$ of 0.025 $\mu$m. This value is very close to the measured experimental displacement (about 0.02 $\mu$m). Note that both the experimental and theoretical central displacements are directed towards the outside of the medium (see white circle arrow in Figures \ref{figElastoPVA}-(A) and \ref{figGreen}-(A)).

To calculate the propagation of the displacements as shear waves, we must consider the transverse dilatation. Indeed, the illuminated zone is a disk of 5 mm in diameter with a thickness of about 40 $\mu m$ (average depth of absorption along Z axis). The dilatation stress is consequently about two orders of magnitude stronger along the transverse direction than along the Z axis: to compute the displacements along time, we neglected the stress along the Z axis and only considered the stress in the transverse direction. We modeled thus the thermoelastic regime in 2D as two opposite forces directed along the Y axis with a depth of 40 $\mu$m and with an amplitude decreasing linearly respectively from 2.5 to 0 mm and from -2.5 to 0 mm. The magnitude of the force along space and time is stored in a matrix, $H_y^{thermo}(y,z,t)$. Displacements along the Z axis are then equal to the convolution between $H_y^{thermo} (y,z,t)$ and $G_{yz}$ \cite{aki1980quantitative}:
\begin{equation}
G_{yz} = \frac{\cos \theta \sin \theta}{4\pi \rho r} (\frac{1}{c_p^2}  \delta_P - \frac{1}{c_s^2} \delta_S +\frac{3}{r^2} \int\limits_{r/c_p}^{r/c_s}{\tau \delta_\tau d\tau})
\label{eq:Gyz}
\end{equation}
where $\delta_P = \delta(t-\frac{r}{c_p})$, $\delta_S = \delta(t-\frac{r}{c_s})$, $\delta_\tau = \delta(t-\tau)$, $c_p$ and $c_s$ are the compression and shear wave speed respectively, $\tau$ is the time, and $\delta$ is a Dirac distribution. The three terms of the equation correspond respectively to the far-field compression wave, the far-field shear wave, and the near-field component.

Using constants $\rho$ = 1000 kg.m$^{-3}$, $c_p$ = 1500 m.s$^{-1}$ and $c_s$ = 5.5 m.s$^{-1}$ (similarly to experimental measurement), displacement maps along the Z axis were calculated 1.0, 1.5, 2.0, 2.5, and 3.0 ms after force application, as illustrated in Figure \ref{figGreen}-(A). The normalized displacement maps present many similarities to the experimental results displayed in Figure \ref{figElastoPVA}-(A), with an initial central displacement directed towards the outside of the medium, and a propagation of three half cycles.