Complex Networks in Phase-Separating Gels. A Computer Simulation Study

1. Introduction

Polymer networks, so-called gels, can undergo two types of phase transitions: volume transitions and phase separation. This study focuses on phase separation, where the gel separates into regions with high and low polymer concentrations. Up to now, research on phase separation in gels has mainly treated surface effects [1,2]. Thus the processes occurring within the bulk of the gel are a welcome object of study for simulations, as we still lack proper methods of physical investigation [1,2]. Comprehensive understanding and prediction of micro phase separation in gels could be of great use for applications where the gel’s microstructure plays a role. Two examples are the development of iongel soft solid electrolytes for solid-state batteries [3,4,5] and the comprehension of the complex deformation behavior exhibited by gels [6].

Motivated by earlier studies [7], which indicated phase separation with a characteristic length scale significantly larger than a model grid size, Peleg et al. [1,2] introduced a simplistic model to mimic the collection of long and short range forces in a gel. The model consists of particles (representing polymer density) interacting via a Lennard-Jones (LJ) potential (short range attraction typifying polymer crystallization) connected through harmonic springs (entropy elasticity) in a two-dimensional simple square lattice structure. Their simulations revealed that, under certain circumstances, the gel phase separates into a high density phase (HDP) and a low density phase (LDP) with the HDP forming a filamentous, foam-like network.

This paper strives to offer further insight into the structure and formation of those HDP networks and to take a step further by generalizing Peleg et al.’s [1,2] model to three dimensions. Of particular interest were the parameters allowing the formation of such structures, their process of formation and their morphology.

2. Materials and Methods

2.1. Model

In the two-dimensional elastic LJ toy model proposed and investigated by Peleg et al. [1,2,8] the gel is represented by particles initially arranged in a two-dimensional square (or hexagonal) lattice. Using LJ units, each particle is attributed a mass of one and is connected permanently to its four initially neighboring particles through weightless harmonic springs. These springs have a spring coefficient, denoted as K, and an equilibrium length of zero. Additionally, pairwise interactions between particles occur via a short range 12-6-LJ potential, which is truncated beyond a cutoff radius $r_c=3\times 2^{1/6}$.

By employing periodic boundary conditions (PBC), the network of harmonic springs functions as long-range repulsive forces, causing the particles to move apart, to their simple cubic lattice positions. Conversely, the LJ forces operate at close range, attempting to "glue" the particles together. For three dimensions this model was expanded to a simple cubic lattice (Figure 1a). Here, each particle is connected permanently to its six initially neighboring particles. The dimension dependent quantities such as distance and temperature were adapted accordingly.

2.2. Molecular Dynamics

The three-dimensional model was studied by a molecular dynamics (MD) simulation in the canonical ensemble. LJ-units are used throughout, i.e., physical units were kept dimensionless. The system was initialized as a simple cubic lattice of N particles and the volume of the simulation box $V=L^3$ was set to give the desired number density $\rho =N/V$. At startup, each particle was attributed a pseudo-random velocity vector $\mathbf{v}$ according to a Maxwell-Boltzmann (MB) distribution $f\left(\mathbf{v}\right)\sim e^{-e_{\mathrm{kin}}/T}$ with specific kinetic energy $e_{\mathrm{kin}}=v^2/2$, while the overall momentum was constrained to be zero. Newton’s equations of motion were integrated by a Velocity Verlet algorithm with an integration time step $\Delta t=0.004$. Temperature was kept constant by periodic velocity rescaling, applied every hundredth time step.

2.3. Force Calculations

The springs were set up in a way, that permanently connected a particle to the nearest periodic-boundary-condition copy of its four/six initial neighbors. In all simulations, spring lengths stayed well below one half of the simulation box length, ensuring consistency. For d dimensions only $d\times N$ springs had to be checked for the spring force calculations. Thus they used only a small part of the total computing time.

The main bottleneck of these large MD simulations were the calculations of the LJ forces. With more than hundred thousand particles and therefore (potentially) several billion LJ force calculations per time step it became imperative to optimize this process. This was achieved by the implementation of neighbor lists. This effectively reduced the number of LJ force calculations per time step by grouping the particles into neighborhoods. The neighborhood of a given particle was made to consist of all other particles, which were likely to interact with it until the next update of the neighbors lists. In the present work, the neighborhood of a particle encompassed all other particles closer than an estimated radius:

$$R_{\mathrm{NL}}=r_c+\Delta R\gtrsim r_c+12\Delta i_{\mathrm{NL}}\Delta t\sqrt{{\displaystyle \frac{2T}{\pi }}},$$

(1)

where $r_c$ denotes the cutoff radius, $\Delta i_{\mathrm{NL}}$ the number of time steps between neighbor list updates, $\Delta t$ the time step, and T the temperature.

These neighbor lists are most effective at low temperatures facing an even particle distribution at low density. At larger temperatures, the neighborhood radius $R_{\mathrm{NL}}$ (1) increases and hence the neighborhoods encompass more particles, which reduce their effectiveness. An uneven particle distribution has a similar effect, increasing the particle count of neighbor lists of all HDP particles drastically. It was chosen to update the list every hundredth time-step $\Delta i_{\mathrm{NL}}=100$. We checked that the criterion (1) ensures that all interactions are properly taken into account for all systems studied (see Appendix A).

As programming language and framework MATLAB was used. All two-dimensional simulations were run and analyzed on a desktop PC (Intel i7-6700 CPU), the three-dimensional simulations were run on the ETH-Zurich Euler cluster. Attempts at using LAMMPS [9] to simulate this model were unsuccessful.

2.4. Methods of Investigation

The evolving HDP networks were of a complex nature. Several methods have been found to describe the gels structure: (i) graphical representation, (ii) minimal particle distribution (MPD), (iii) spring length distribution (SLD), and (iv) radial shell density distribution (RSD).

(i) Figure 1b and c show a graphical representation of the trajectories of the simulated gel. The connecting springs were drawn with a limited number of points per line. This had the effect of rendering extended springs transparent and contracted springs opaque. In two dimensions, springs smaller than $2\times 2^{1/6}$ were colored orange highlighting the HDP. In three dimensions, foam-like structures consisting of faces and filaments evolved. To increase the visibility of these structures, each spring was colored according to its spacial orientation. Over all this had the effect of highlighting faces and coloring them according to their orientation.

(ii) To determine the MPD of a gel state, for each particle the distance to its closest neighbor was measured. Then a binning algorithm was applied, with thousand equally spaced binning intervals ranging from zero to one fifth of the simulation box length.

(iii) To construct the SLD, this same binning algorithm was used on all current spring lengths instead of the minimal radial particle distances.

(iv) The RSD, which is the radial distribution function not normalized by the bulk density, is usually calculated by measuring all pairwise particle distances, then adding, binning and normalizing them. With a gel of $N={60}^3$ particles this would have meant calculating and binning ${60}^6$ distances. To keep the calculation effort reasonable, for each particle one percent of the other particles were chosen at random to then measure the distances. The binning was carried out with 2500 binning intervals over half of the simulation box length and therefore with the same binning interval length, or more precisely the same thickness of the spherical binning shell $dr$ as for the MPD and SLD. The resulting distribution was then normalized by the volume of the binning shell ($4\pi r^{2}dr$) and by the number of sample distances per particle ($0.01\times N$).

To monitor the evolution of the structure over time, several measures were taken to express the gels structure in one number per time frame. This process is sketched in Figure 2. Any two particles closer to each other than a critical distance of $d_{cd}=1.1\times 2^{1/6}$ were declared HDP particles (see Figure 2). The simulation box was then regularly separated into small cubes of edge length $s_b=2$. To cancel out random agglomerates, such a small cube needed to encase more than four HDP particles to be declared full. Thereafter a new grid of big cubes with edge length $s_B=2s_b$ was created - each big cube therefore consisting of eight small cubes in three dimensions or four in two dimensions. A big cube was now declared to be a surface cube, if at least one, but not all of its small cubes were full. The above values for agglomerate threshold, the size of a small box, the number of small cubes making up a big cube and the upper and lower limit for a big cube to be a surface cube were adjusted graphically to fit the HDP structures visible by eye as good as possible.

By this algorithm the following three numbers were obtained: the number of HDP particles $N_{\mathrm{HDP}}$, the number of full (small) cubes $b_{\mathrm{full}}$, and the number of (big) surface cubes $B_{\mathrm{surf}}$. From these, five quantities of interest could be obtained for two and three dimensions d: namely the relative HDP surface area

$$a_{\mathrm{HDP}}={\displaystyle \frac{A_{\mathrm{HDP}}}{L^{d-1}}}={\displaystyle \frac{1}{L^{d-1}}}\times d{B}_{\mathrm{surf}}\times s_B^{d-1},$$

(2)

the fraction of HDP particles (only the particles inside of the small cubes were counted)

$$n_{\mathrm{HDP}}={\displaystyle \frac{N_{\mathrm{HDP}}}{N}},$$

(3)

the HDP volume fraction

$$v_{\mathrm{HDP}}={\displaystyle \frac{V_{\mathrm{HDP}}}{L^d}}={\displaystyle \frac{1}{L^d}}\times b_{\mathrm{full}}\times s_b^d,$$

(4)

the density of the HDP

$${\rho }_{\mathrm{HDP}}={\displaystyle \frac{N_{\mathrm{HDP}}}{V_{\mathrm{HDP}}}}={\displaystyle \frac{N_{\mathrm{HDP}}}{b_{\mathrm{full}}\times s_b^d}},$$

(5)

and lastly the thickness of the HDP filaments

$${\delta }_{\mathrm{HDP}}={\displaystyle \frac{(d-1)V_{\mathrm{HDP}}}{A_{\mathrm{HDP}}}}.$$

(6)

This Eq. (6) approximates the filament as an infinitely long cylinder, where the total surface is dominated by the curved surface area. There the ratio of volume to surface is

$${\displaystyle \frac{V_{\mathrm{cy}}}{A_{\mathrm{cy}}}}={\displaystyle \frac{\pi r^{2}h}{2\pi rh}}={\displaystyle \frac{r}{2}}\mathrm{leading}\mathrm{to}r={\displaystyle \frac{2V_{\mathrm{cy}}}{A_{\mathrm{cy}}}}$$

(7)

for the radius r. This is the ratio used in Eq. (6) to get a rough measure of thickness.

The same argument can be made in two dimensions, taking a rectangle instead of a cylinder.

In addition to these structural measures, the total static energy of the system $E_{\mathrm{tot}}=E_{LJ}+E_{\mathrm{K}}$ was calculated as the sum of the total pairwise LJ energy

$$E_{LJ}=4\sum _{i=1}^N\sum _{j>i}^N\left({\displaystyle \frac{1}{r_{ij}^{12}}}-{\displaystyle \frac{1}{r_{ij}^6}}\right)$$

(8)

and the total spring energy

$$E_K={\displaystyle \frac{K}{2}}\sum _i^N\sum _{j:j⌣i}^d{r}_{ij}^2.$$

(9)

Here $r_{ij}$ denotes the scalar distance between particles i and j, and $i⌣j$ symbolizes a spring connection between them.

3. Results

3.1. Evolution of a HDP Network

Under certain conditions, the gel could be observed to phase separate into a complex network of HDP, enclosing regions of LDP. Figure 3 shows this at the example of a two- and a three-dimensional gel. In both examples, small HDP precipitates could be seen from the start. Over time they grew and connected to each other to form a complex HDP network.

To monitor this evolution, the three-dimensional gel’s trajectories were analyzed according to the measures established in Section 2.4. Figure 4a shows particle fraction $n_{\mathrm{HDP}}$, volume fraction $v_{\mathrm{HDP}}$ and density of the HDP. All three measures showed a rapid increase at the beginning which slowed down over time until asymptotic behavior was prevalent. The final volume occupied by the HDP was little less than $3\%$ concentrating $73\%$ off all particles to a density of ${\rho }_{\mathrm{HDP}}=0.81$ inside the HDP.

The evolution of the relative HDP surface area $a_{\mathrm{HDP}}$ is shown in Figure 4b. Two stages can be observed. First, as HDP precipitated and precipitates grew larger than the agglomerate threshold, the relative HDP surface area increased. Later, after reaching a maximum of $a_{\mathrm{HDP}}=1.44$ (HDP surface area was $44\%$ larger than one $L^2$ side area of the surface cube) between $t=500$ and $t=1000$, the relative HDP surface area decreased again, until transitioning to asymptotic behavior towards end of the simulation. Relative HDP surface area $a_{\mathrm{HDP}}$ and HDP volume fraction $v_{\mathrm{HDP}}$ are linked over Eq. (6) to give the HDP filament thickness ${\delta }_{\mathrm{HDP}}$. This HDP filament thickness decreased to a minimum just until the onset of HDP network formation visible by eye (see Figure 3) and then increase to an asymptotic value of ${\delta }_{\mathrm{HDP}}=5.5$ (simulation box length was about $L=130$).

The evolution of the gels static energies showed a decrease in total LJ energy $E_{\mathrm{LJ}}$ with a lesser increase in total spring energy $E_{\mathrm{K}}$ leading to a decrease in total static energy $E_{\mathrm{tot}}$ over the course of the simulation (see Figure 4c). After an erratic first phase, the changes in all energies slowed down to become asymptotic towards the end of the simulation.

To further analyze the morphology of the evolving structure, MPD, SLD and RSD were consulted (see Figure 5). All three distributions showed a peak around the LJ minimum indicated by a black dashed-dotted line at a distance (or length or radius resp.) of $2^{1/6}$. The dashed blue line in the MPD separates HDP from LDP particles. Looking closely at the tip of the MPD at a distance of $d_{\mathrm{MPD}}=2^{1/6}$, a shift to lower distances can be noticed, indicating an increasing HDP density. After $t=500$ two peaks at $d_{\mathrm{MPD}}=1.9$ and around $d_{\mathrm{MPD}}=4$ became distinct - the right-hand-side peak moving towards larger distances, which indicating LDP regions of decreasing density. In the SLD, peaks also grew more distinct over the course of the simulation. At all times, the right-hand-side of the SLD showed a steep decrease for the number of large spring lengths $s_{\mathrm{SLD}}>5$. In this run SLD at $s_{\mathrm{SLD}}=2^{1/6}$ could be observed to decrease over time from $5\times {10}^4$ at the beginning to $2\times {10}^4$ after $t=4000$. As more and more particles joined the HDP over time, the radial shell density grew, and its peaks became more distinct. Drawing the distances of the second to fourth closest neighbors of a bcc and fcc structure with minimal distance $2^{1/6}$ into the RSD, suggests that a large fraction of the HDP is fcc.

3.2. Parameter Variations

Originating from the parameter set of the two- and three-dimensional examples treated in the previous Section 3.1, the parameter space spanned by density $\rho$, temperature T and spring coefficient K was further explored. The influence of $\rho$, T and K could be investigated separately by choosing by choosing one of these parameter to vary, while the other two were kept constant.

The two-dimensional (2D) gel model simulations were run with $N={100}^2$ particles originating from ${\rho }_{2\mathrm{D}}=T_{2\mathrm{D}}=K_{2\mathrm{D}}=0.3$, the three-dimensional (3D) simulations with $N={60}^3$ particles, originating from ${\rho }_{3\mathrm{D}}=0.1$ and $T_{3\mathrm{D}}=K_{3\mathrm{D}}=0.3$. Figure 6 shows the graphical representation of these simulations after a time of $t=4,000$.

The influence of varying the density on the morphology of the end state of the simulations is shown in Figure 6a. At lowest densities, no phase separation was observed, HDP did not precipitate. At slightly higher densities ${\rho }_{2\mathrm{D}}=0.2$ and ${\rho }_{3\mathrm{D}}=0.05$, large thin structures could be observed: thin HDP filaments with large LDP holes in-between in the 2D simulations and large thin faces, normal to the initial simple cubic spring lattice, in the 3D simulations. Increasing the density further had the effect of making the HDP network bulkier, especially at the intersections of the filaments (2D simulations) or the faces (3D simulations) and of lowering the observed length scale of the network. At the highest densities ${\rho }_{2\mathrm{D}}=0.6$ and ${\rho }_{3\mathrm{D}}=0.4$, HDP made up the larger part of the simulation box, with only few LDP holes left.

At too low temperatures $T_{2\mathrm{D}}=T_{3\mathrm{D}}\le 0.1$, HDP precipitated but the small precipitates did not connect to a larger network (see Figure 6b). Increasing the temperature to $T_{2\mathrm{D}}=0.4$ and $T_{3\mathrm{D}}=0.5$, lead to an increase in perceived length scale of the network – fewer, but larger faces in the 3D simulation and larger thicker filaments in the 2D one. Especially in the 2D simulation, the filaments appeared less compact, fuzzier than in the ’origin’ simulations at $T_{2\mathrm{D}}=T_{3\mathrm{D}}=0.3$. For larger temperatures $T_{2\mathrm{D}}\ge 0.5$ and $T_{3\mathrm{D}}\ge 1$, no phase separation was observed.

If the spring coefficient K is zero, the springs are inactive and the gel becomes a LJ fluid instead, which is already well researched [10,11]. Therefore too low values of $K<0.05$ were not of interest. At the lowest investigated spring coefficient $K_{2\mathrm{D}}=0.05$ of the 2D model, thin filaments could be observed, connecting junks of HDP (see Figure 6c). The color of the observed HDP was only faintly orange, indicating the presence of many springs longer than $2\times 2^{1/6}$. The 3D gel at $K_{3\mathrm{D}}=0.1$ showed only HDP precipitates with no apparent connections.

Moving to larger spring coefficients $K_{2\mathrm{D}}=0.1$, $0.3$, the orange coloration of the HDP in the 2D gels got stronger, indicating more short springs inside the HDP. Additionally, the HDP filaments grew thicker, with now only slightly bulkier intersections. A similar morphology could be observed in the 3D gel at $K_{3\mathrm{D}}=0.3$: a network of faces connected at slightly bulkier intersections.

At large spring coefficients $K_{2\mathrm{D}}=0.5$, $0.9$ the filaments of the 2D gel became very straight, the LDP space in-between changing in shape from roundish to edgy. Less but larger LDP areas could be seen. The LDP bubbles of the 3D gel behaved in a similar fashion: larger, but still slightly curved faces at $K_{3\mathrm{D}}=0.5$ turning into large planar constricted faces at $K_{3\mathrm{D}}=0.9$. For $K_{2\mathrm{D}}\ge 1$ and $K_{3\mathrm{D}}\ge 1$ no phase separation was observed.

3.3. Phase Diagrams

To create an overview, at which combinations of density $\rho$, temperature T and spring coefficient K phase separation and network formation could be found in the three-dimensional model, many small simulations of $N={30}^3$ and $N={40}^3$ particles were started, and run until $t=4000$ (${10}^6$ time steps), or until HDP network growth could be seen. As simulation times in three dimensions were quite high (a few weeks at this size) these simulations were limited to a plane of constant density $\rho =0.1$.

Figure 7 shows the phase diagram of the three-dimensional model, indicating where phase separation occurred (pink dots/regions I and II), if the HDP formed a network (blue circles/region II) or if neither thing happened (light blue stars/region III). The classification of the gel simulations into these three regions was carried out by the author, judging the graphical representation by eye.

The same approach was used to map the HDP network formation of the two-dimensional model. Simulations of $N={100}^2$ and $N={150}^2$ particles were used to investigate the gel on three planes of constant density $\rho =0.1$, $0.2$ and $0.3$. Figure 8 shows the resulting three phase diagrams.

For all three planes of density, a similar shape of the three regions was found, the HDP network region (II) more askew to the left, in the lower density phase diagrams. Lowering the density had a larger impact on the upper spring coefficient limit, than on the temperature limit of region II. The phase separation regions (I and II) of the two-dimensional model phase diagram with $\rho =0.1$ was smaller than phase separation regions of its corresponding three-dimensional model phase diagram (Figure 7).

4. Discussion

4.1. Static Energy Balance between HDP and LDP

The phase separation of the gel can be explained with the competing LJ and harmonic energies. To do so, the toy model was simplified to a one dimensional string of N LJ particles connected through harmonic springs with one group of $n+1$ HDP particles and the rest LDP particles (see Figure 9). In this simplification, the effect of thermal vibrations is neglected - all HDP particles are assumed to be precisely in the LJ minimum at $r_{\min }=2^{1/6}$ and $E_{\min }=-1$. Therefore all HDP springs are of equal length $s_1=s_2=...=s_n=r_{\min }$. Further we assume the LDP springs to be spaced equally and to be larger than the LJ cutoff radius. Thus $s_{n+1}=s_{n+2}=...=(L-n{s}_1)/(N-n)>r_c$.

The LJ energy of the n HDP particles then is $E_{\mathrm{LJ}}=\left(-1\right)\times n$ as only the LJ energies of the HDP particles is nonzero. The energy of the harmonic springs is

$$\begin{array}{cc}\hfill E_{\mathrm{hrm}}& ={\displaystyle \frac{K}{2}}\sum _{i=1}^n{s}_i^2+{\displaystyle \frac{K}{2}}\sum _{j=n+1}^N{s}_j^2={\displaystyle \frac{K}{2}}\left(n{s}_1^2+(N-n){\displaystyle \frac{{(L-n{s}_1)}^2}{{(N-n)}^2}}\right)\hfill \\ \hfill & ={\displaystyle \frac{K}{2}}\left(n{s}_1^2+{\displaystyle \frac{{(L-n{s}_1)}^2}{(N-n)}}\right).\hfill \end{array}$$

(10)

The graph of these two energies and their sum $E_{\mathrm{tot}}=E_{\mathrm{LJ}}+E_{\mathrm{K}}$ is given in Figure 9b. For $N=100$ particles, a spring coefficient of $K=1$ a density of $\rho =0.75$ and a box length of $L=N/\rho =133$ the total static energy has a minimum at $n=85$.

Varying density and spring coefficient, a surface of minimal total static energy was found - giving for each combination of $\rho$ and K the optimal HDP particle fraction $n_{\mathrm{HDP}}$ (see Figure 9c). This shows an energetic motivation for phase separation, driving the system towards an equilibrium state of minimal energy.

Considering one alteration of the model, the observed growth of the network can be explained: If we imagine not one, but several groups of HDP particles, for each new group created, one LJ bond is lost and replaced with a stretched spring, raising the total energy. The observed surface optimization is the reverse of this imagined process - decreasing the total energy through replacing long springs with LJ bonds.

4.2. Formation and Dissolution of HDP

To explain the influence of the temperature, and why sometimes phase separation occurred without the formation of a HDP network, let us consider an other one-dimensional model (see Figure 10a): One, mobile LJ particle is connected with two springs of lengths $s_1$ and $s_2$ to its neighbors which are at a fixed distance of $2s_0=s_1+s_2$ to each other. A short calculation reveals the spring energy to be

$$E_{\mathrm{K}}={\displaystyle \frac{K}{2}}(s_1^2+s_2^2)={\displaystyle \frac{K}{2}}{(r^2+(2s_0-r))}^2)=K(r^2-2r{s}_0+2s_0),$$

(11)

where r denotes the distance of the moving particle from its left-hand-side neighbor. Together with the LJ potential, it follows that for $r<s_0$, the total potential of the mobile particle is

$$E_{\mathrm{tot}}=K(r^2-2r{s}_0+2s_0)+4\left({\displaystyle \frac{1}{r^{12}}}-{\displaystyle \frac{1}{r^6}}\right).$$

(12)

Both energies $E_{\mathrm{K}}$ and $E_{\mathrm{LJ}}$ and the resulting total energy $E_{\mathrm{tot}}$ are plotted in Figure 10b.

Figure 10. (a) Sketch of a LJ particle connected to two other fixed particles over harmonic springs with equilibrium length zero. The mobile particle is at a distance of $r<s_0$ to its left neighbor. (b) Harmonic spring energy $E_{\mathrm{K}}$, LJ energy $E_{\mathrm{LJ}}$ and total energy $E_{\mathrm{tot}}$ of the mobile particle for $r<s_0$ with $K=0.25$ and $s_0=5$.

Figure 10. (a) Sketch of a LJ particle connected to two other fixed particles over harmonic springs with equilibrium length zero. The mobile particle is at a distance of $r<s_0$ to its left neighbor. (b) Harmonic spring energy $E_{\mathrm{K}}$, LJ energy $E_{\mathrm{LJ}}$ and total energy $E_{\mathrm{tot}}$ of the mobile particle for $r<s_0$ with $K=0.25$ and $s_0=5$.

It can be observed, that in order to get to the LJ minimum, which is only a local minimum in the regarded example, the mobile particle needs to overcome a considerable energy barrier.

If we would extend this model to the right-hand-side, considering a second a third or more mobile particles, those particles could sometimes aid this first one with reaching the left LJ minimum by approaching or by accelerating it to the left, and other times keep it from reaching the LJ minimum by doing the opposite.

Exact calculations would need to be conducted through a statistical approach, considering the possible vibrations of the gel. For now the intuition must suffice: It is energetically adverse for one particle to reach the other if it is moving alone. It needs more energy, through either a higher temperature, or through neighboring particles vibrating collectively, which can be achieved by slightly increasing the spring coefficient. If there is too much thermal energy or if those neighbors vibrate too much, the particle is in danger of being ripped out of the LJ bond again.

This intuition can be applied to the phase diagrams found for the two- and three-dimensional gel models (see Figure 7 and Figure 8). At low temperatures and low spring coefficients, particles didn’t have enough vibrational energy to reach each other and form a network. Only small precipitates form (Region I). This could be remedied by reinforcing the connection between the particles by increasing the spring coefficient, to make them more likely to vibrate collectively and thus enabling the precipitates to form a network (region II). Too large spring coefficients then had the reverse effect, forcing the particles apart and dissolving the HDP again (region III). Another possibility of connecting the precipitates from region I was to increase the temperature directly, giving the particles more energy to bridge the gap between precipitates. Too large temperatures then dissolved the HDP again.

4.3. Future Works

Shorter simulation times, allowing more or larger simulations, would be very welcome for future projects. The most promising alternative to the rather slow MATLAB code was found in rewriting the simulation in Julia [12]. This way, a simulation speed improvement of factor nine could be achieved. More might be possible by improving the neighbors lists or by parallelizing the calculations.

Future works treating the same model could quantify the influence of simulation parameters $\rho$ ,T and K on morphology using the final measures from Section 2.4 (HDP surface fraction $a_{\mathrm{HDP}}$, HDP thickness ${\delta }_{\mathrm{HDP}}$, ...), to move away from describing the observed change in morphology from images as done in Section 3.2. This would remove subjectivity and allow possible automatization to e.g. find the parameters for the largest surface area.

For the same reason, fix, non-subjective criteria should be developed for associating gels to regions I, II and III.

One possibility would be, to use the small boxes covering the HDP volume (established in Section 2.4). Region II could be defined as a percolating structure of small boxes. This could be implemented with a forest fire algorithm like in [13,14]. Other suitable criteria could be found, taking a certain minimal HDP particle fraction combined with a minimal HDP thickness. Either would enable an automated mapping of phase diagrams like Figure 7 and Figure 8.

Instead of running simulations with constant parameters, outgoing from an initial state, simulations could be started from the final/asymptotic state of a simulation of other parameters. This would serve to further investigate the temperature hysteresis Peleg et al. [1] found, to see if a hysteresis could also be observed in the other two parameters ($\rho$ and K), and if the same phase diagram would develop eventually.

The prediction for the optimal HDP particle fraction of a one-dimensional gel model (Section 4.1) could be compared to simulations of a one-dimensional gel model, and generalized to also include the influence of temperature.

To get closer to reality, the model could be altered by e.g. using a distribution of spring coefficients instead of a uniform K or by randomizing the structure through e.g. cutting/disabling some of the springs of the initial simple cubic lattice, as in [8]. Another approach would be to insert a second species of solvent particles.

The two-dimensional model of Peleg at al. [1,2] has already been compared to experimental results from investigating biopolymer mixtures in confined geometries [15] and found to be "in good agreement" with their findings.

An other direct link to observable gels might be achieved through changing the three-dimensional model to simulate a surface, by changing one dimension from periodic to non-periodic with a boundary against vacuum or LJ fluid. The resulting surface could then be compared to scanning electron microscopy (SEM) images or atomic force microscopy (AFM) images of gel surfaces like they were obtained in [16,17,18,19].

Alternatively, the resulting structure of larger three-dimensional simulations could be compared to gels investigated with micro-computed tomography (micro-CT) [20].