## 3 Physics-Based Particle Simulation Approach  
### Overview of the DPD Simulation Method  
DPD was employed as the primary simulation technique to model the mesoscale behavior of emulsifier-oil-water systems. DPD is a coarse-grained molecular dynamics approach that enables the simulation of larger spatial and temporal scales by grouping atoms into beads and applying simplified interaction rules. Each DPD particle represents a cluster of molecules, and their interactions are governed by Newtonian dynamics incorporating conservative, dissipative, random, and bonding forces [36]. This method is particularly suited for capturing the hydrodynamic behavior of multiphase systems such as emulsion.  
Figure 1: Simplified methodology workflow for multiscale simulation of HIPE formation.  
In DPD simulations, the total force acting on a particle arises from the interplay of several components, primarily \(F_{ij}^{C}\), the dissipative force \(F_{ij}^{D}\), and the random force \(F_{ij}^{B}\). In addition, a spring force \(f_{j}^{S}\) was also taken into consideration, as described below [37]:  
\[f=\sum_{i\
eq j}\left(F_{ij}^{C}+F_{ij}^{D}+F_{ij}^{R}\right)+f_{i}^{S} \tag{1}\]  
where \(F_{ij}^{C}\) is the conservative force (repulsion force) and models excluded volume effects, \(F_{ij}^{D}\) is the dissipative force, introducing friction to maintain local momentum conservation, \(F_{ij}^{R}\) is the random force, representing thermal fluctuations, and \(f_{j}^{S}\) is the spring force which accounts for bonding interactions between connected beads. The conservative force describes the maximum repulsive interaction between beads \(i\) and \(j\), given as [38]:  
\[F_{ij}^{C}=\begin{cases}a_{ij}(1-r_{ij})r_{ij}&\text{if }|r_{ij}|<1\\ 0&\text{if }|r_{ij}|\geqslant 1\end{cases} \tag{2}\]  
where \(a_{ij}\) is the maximum repulsion parameter derived from the solubility parameters of the beads (typically calculated via atomistic MD). \(r_{ij}\) is the distance between the beads, and \(r_{ij}\) is the unit vector along the bead-bead separation. The dissipative force introduces friction proportional to the relative velocity between beads [38]:  
\[F_{ij}^{D}=\begin{cases}-\gamma\omega^{D}r_{ij}(\hat{r_{ij}}v_{ij})r_{ij}& \text{if }|r_{ij}|<1\\ 0&\text{if }|r_{ij}|\geqslant 1\end{cases} \tag{3}\]  
where \(\gamma\) is a friction coefficient, \(v_{ij}\) is the relative velocity of beads \(i\) and \(j\), and \(\omega^{D}r_{ij}\) is a weight function dependent on the inter-particle distance. The random force accounts for thermal noise and is expressed as [38]:  
\[F_{ij}^{R}=\begin{cases}\sigma\omega^{R}r_{ij}\theta_{ij}r_{ij}&\text{if }|r_{ij}|<1\\ 0&\text{if }|r_{ij}|\geqslant 1\end{cases} \tag{4}\]  
where \(\theta\) is a random flocculation variable between 0-1, \(\sigma\) is the amplitude of the noise, and \(\omega^{R}r_{ij}\) is a weighting function related to \(\omega^{D}r_{ij}\)[38]:  
\[\omega^{D}(\mathbf{r})=\left[\omega^{R}(\mathbf{r})\right]^{2}=\begin{cases}(1 -\mathbf{r})^{2},&\left|r_{ij}\right|<1\\ 0,&\left|r_{ij}\right|\geqslant 1\end{cases} \tag{5}\]  
The amplitude of the random force (\(\sigma\)) and the dissipative force coefficient (\(\gamma\)) are related to the system's temperature (\(T\)) through the fluctuation-dissipation theorem:  
\[\sigma^{2}(r)=2\gamma k_{B}T \tag{6}\]  
where \(k_{B}\) is the Boltzmann constant while \(T\) is the absolute temperature.  
All simulations were performed using dimensionless DPD units with a reduced cutoff radius \(Rc=1\), mapped to 3.23 A. The simulation timestep was 0.005\(tc\), and the system temperature was maintained at 65 \({}^{\circ}\)C in the NVT ensemble over 1,000,000 steps [39].