Analysis of bijel formation dynamics during solvent transfer-induced phase separation using phase-field simulations

Abstract

Solvent transfer-induced phase separation (STrIPS) is a versatile and scalable method for the fabrication of bicontinuous interfacially jammed emulsion gels (bijels), a soft material with ample potential in membrane separation and energy storage. STrIPS, however, does not facilitate in situ observation of bijel formation, rendering it difficult to relate the phase separation process to the resulting bijel morphology. To address this issue, a phase-field (PF) model is presented that captures the principle components of STrIPS through relatively simple methods. Comparison with experimental results demonstrates that the PF model accurately describes STrIPS, reproducing morphological trends in a complex experimental system. Simultaneously, the model provides insight into the origin of such morphological features, relating the formation of the bijel structure in STrIPS to relatively slow dynamics of phase separation compared to solvent diffusion. Visualising phase separation during the entire STrIPS process, the PF model can therefore help elucidate the formation mechanisms underlying current STrIPS bijel morphologies, in addition to aiding the fabrication of novel ones.

---

1 Introduction

The stabilisation of emulsion droplets with solid particles has been a fixture in colloid science over the past few decades.^1–4 More recently, this concept was extended to include bicontinuous emulsions, giving rise to the so-called bicontinuous interfacially jammed emulsion gel (bijel).^5,6 Owing to its unique morphology, the bijel has ample potential for application in a wide variety of industrially relevant fields, such as membrane separation,^7,8 coatings,⁹ energy storage^10,11 and catalysis.^12,13

While multiple methods for the fabrication of bijels exist,^6,14,15 solvent transfer-induced phase separation¹⁶ (STrIPS) offers unique advantages in terms of composition versatility^17,18 and process scalability.¹⁹ During STrIPS, the removal of a solvent induces phase separation of a liquid precursor mixture additionally containing water, oil and interfacially active nanoparticles. The nanoparticles subsequently attach to the interface between the newly formed immiscible phases, eventually creating a jammed monolayer that kinetically arrests the structure and prevents further coarsening.

By employing a solvent to induce phase separation, the range of suitable precursor liquids, previously limited to only a few binary liquid pairs with select thermal properties,^6,14,20,21 is greatly expanded. Simultaneously, the inclusion of the solvent enables the transition from a purely batch-based to a continuous synthesis process.

While STrIPS has these distinct advantages over other fabrication methods, it does not allow for facile in situ observation of bijel formation through conventional techniques (e.g. confocal microscopy). As such, current research often attempts to establish relationships between the initial state of the precursor mixture and the final bijel morphology, while the intermediate processes remain unclear. Recent work has partially addressed this issue through time-controlled solidification of a polymerisable oil,^7,8 but this remains limited in scope due to the restrictive selection of suitable precursor components.

Instead, simulations can provide useful insights into bijel formation during the STrIPS process. Relating the properties of the precursor mixture's constituent components to the observed phase separation dynamics can help rationalise previously obtained experimental results, in addition to aiding the realisation of new, specifically tailored bijel morphologies.

Simulations have always been prevalent in research surrounding the bijel, being first discovered through numerical simulations by Stratford et al.⁵ prior to its experimental realisation by Herzig et al. two years later.⁶ Later work explored the influence of different liquid volume fractions,²² particle properties such as arbitrary size, shape or wettability^23,24 and substrate wetting effects.²⁵ While these simulations provided invaluable insights into the fundamentals of bijel formation, they are generally based on either the Lattice-Boltzmann (LB)^{5,22,24,26–28} or dissipative particle dynamics (DPD)²⁹ methods, severely limiting the spatial and temporal scales over which simulations can be performed (micrometers over the course of milliseconds). For STrIPS these limitations are especially problematic, as the length and time scales of interest are set by the macroscopic molecular diffusion of the solvent (hundreds of micrometers over the course of seconds).

To simulate the formation of bijels on these larger scales, phase-field (PF) models are a promising alternative. Rather than the discrete LB or DPD methods, PF models work by solving partial differential equations of continuous fields that represent the different components of the system, with interactions that follow from mean-field approximations. In this manner, PF models can access time and length scales that are unfeasible for either LB or DPD. PF models have already found extensive employment in the simulation of phase separation processes, particularly related to the formation of polymeric membranes through phase inversion.^30–34 Moreover, a combined PF/Brownian dynamics model has been used to investigate bijel formation in thin-film geometries.^23,25,35 The explicit inclusion of particles and their interactions, however, still renders these types of simulations impractical for the large scales associated with the STrIPS process.

In this work, a PF model is introduced that captures the key characteristics of the STrIPS process while remaining computationally inexpensive. To this end, two distinct phase-field variables are evolved that together capture the physics of the oil/water/solvent system. In this model, the arrest of phase separation does not come about through the modelling of explicit particles. Rather, a local measure of the interfacial tension is continuously calculated and phase separation arrested when it exceeds a certain threshold value. In this way, the process of particles attaching to and jamming the formed interfaces is effectively replicated.

After the description of the underlying theory and subsequent model verification, the influence of different model parameters on the structure of the STrIPS bijel is investigated. The parameters controlling the relative dynamics of phase separation and solvent diffusion are found to be especially crucial and their role in the formation process is further elucidated. Finally, model results are compared with experimentally obtained STrIPS bijel morphologies, showing strong qualitative agreement. In particular, the potential of the PF model is demonstrated by readily capturing trends in the morphological features of the STrIPS bijel based on different experimental conditions.

2 Theory

The STrIPS process can be separated in three distinct components: solvent diffusion, phase separation and interfacial nanoparticle jamming. While the former two can be described via the same set of equations,³⁰ it is simpler to treat them as separate processes. In this work, this is achieved by simultaneously evolving two different phase-fields, one which represents the oil/water mixture and one that represents the solvent. By making the oil/water interaction dependent on the value of the solvent field, the physics of the oil/water/solvent system can be replicated.

2.1 Oil/water phase separation

In this work, the phase separation of the oil/water mixture is described in accordance with Cahn–Hilliard theory.³⁶ Herein, the thermodynamic properties of the mixture can be derived from its free energy functional F, which in addition to the bulk free energy density f also contains gradient energy contributionswhere V is the system volume, ϕ is chosen to represent the volume fraction of the oil phase and κ is the gradient penalty coefficient. From the free energy functional, the chemical potential of the oil phase can subsequently be defined as the variational derivative.Assuming mass transport occurs solely through diffusion, the chemical potential of eqn (2) can be combined with a general continuity expression to yield the Cahn–Hilliard equation aswhere for the latter equality it is assumed that the mobility of the oil phase M remains constant.

In this work, the oil/water mixture is treated as a simple regular solution with a free energy density given by

f = α(ϕ ln ϕ + (1 − ϕ)ln(1 − ϕ) + χϕ(1 − ϕ)) (4)

where α = k_bT/V is the scale of the free energy density with the Boltzmann constant k_b and the absolute temperature T. Here, 1 − ϕ thus represents the volume fraction of the water phase and χ is the oil/water interaction parameter. By coupling this interaction parameter to the solvent phase-field, solvent transfer-induced phase separation of the oil/water mixture can be achieved.

2.2 Solvent diffusion

Separated from the oil/water mixture, the phase-field of the solvent is evolved via Fick's second law of diffusionwhere ϕ_s and D_s represent the solvent's volume fraction and diffusion constant, respectively. The coupling of the solvent field to the oil/water field is achieved by making the oil/water interaction parameter χ linearly dependent on the value of the solvent volume fraction ϕ_swhere χ₀ is the oil/water interaction parameter without any solvent present and χ_c is the critical value of the interaction parameter above which phase separation occurs. For the regular solution description of the free energy density, this entails χ_c = 2. Additionally, the solvent volume fraction ϕ_s is scaled with respect to ϕ_sc, which is a chosen value of the solvent field where χ = χ_c and phase separation thus starts to occur. By scaling in this way, the chosen value of ϕ_sc effectively becomes a measure of the strength of the solvent. This is illustrated in Fig. 1, showing the phase diagrams of the oil/water mixture with respect to the solvent field for different values of ϕ_sc.

Fig. 1 Phase diagrams of the oil/water mixture with respect to the solvent volume fraction ϕ_s for different values of ϕ_sc. The solid and dashed lines are the binodal and spinodal curves, respectively, calculated using eqn (4) and (6) for χ₀ = 3.0.

In these phase diagrams, the two-phase region is located under the calculated binodal curve (solid line). Choosing lower values of the critical solvent fraction ϕ_sc shrinks this region and shifts it down in the phase diagram. As such, the oil/water mixture stays homogeneous at lower solvent values, representing an increase in the solvent strength.

Closer observation of the phase diagram in Fig. 1 reveals one of the main considerations of modelling the diffusion of the solvent separately from the phase separation of the oil/water mixture. Namely, as these two processes are solely coupled through the interaction parameter χ(ϕ_s), there is no condition that enforces that ϕ + ϕ_s = 1. Consequently, this kind of coupling renders it difficult to physically interpret ϕ_s and ϕ as true fractions within the same system.

Rather, as a means to control the interaction parameter between the phase-separating species, the solvent fraction ϕ_s can be considered the STrIPS analogue of the absolute temperature T for the temperature induced phase separation of binary polymer mixtures.^31,32 An alternative, perhaps more physical, interpretation is to treat ϕ_s as the solvent fraction of a composite oil/solvent component ϕ. While this interpretation was originally introduced for the PF modelling of non-solvent induced phase separation of ternary polymer solutions,³⁴ a slightly different formulation of the relation between χ and ϕ_s in eqn (6) would also allow it to be applied here.

2.3 Nanoparticle interfacial jamming

In STrIPS, phase separation is eventually arrested due to the interfacial jamming of adsorbed nanoparticles. However, explicitly modelling nanoparticles is complex and computationally expensive, limiting the scale over which simulations can be performed. Instead, the influence of the nanoparticles can be implicitly included by taking into account the driving force of their interfacial adsorption: the interfacial tension σ. In STrIPS the diffusion of solvent and subsequent formation of interfaces due to phase separation is associated with an increase in the interfacial tension. If the interfacial tension exceeds a certain threshold σ_c, particles are irreversibly adsorbed to the interface and jamming occurs. For a planar 1D interface, the interfacial tension can be calculated in accordance with³⁷For a 2D system this can be extended by defining a local measure of the surface tension as the geometric mean of the surface tensions calculated in perpendicular directionshere, the interfacial jamming of nanoparticles is thus replicated by calculating the local interfacial tension σ in accordance with eqn (8) and freezing the phase-field of the oil/water mixture when it exceeds a critical value σ_c. The details regarding the calculation of σ via eqn (8), in addition to its relation with the solvent fraction ϕ_s, can be found in Section S2 of the ESI.†

3 Methods

3.1 Nondimensionalisation

In the PF model, the governing equations of phase separation and solvent diffusion, respectively eqn (3) and (5), are numerically solved in a dimensionless form. To this end, the following dimensionless variables are introducedwhere l is a characteristic length present in the simulation and D_s the diffusion coefficient of the solvent. Thus, the time is respectively scaled to the characteristic diffusion time of the solvent over this length, τ = l²/D_s, such that eqn (5) becomesScaling with respect to the solvent allows the PF model to capture the dynamics of both the solvent and the oil/water components, while simultaneously concentrating the remaining nondimensional parameters in the Cahn–Hilliard equationwhere Ψ = f/α is the dimensionless version of the free energy density in eqn (4). In this form, M̂ = Mα/D_s is the relative mobility of the oil/water mixture with respect to the solvent, whereas [small kappa, Greek, circumflex] = κ/(l²α) sets the relative energy penalty of creating compositional gradients, compared to the bulk free energy.

3.2 Numerical implementation

The nondimensionalised eqn (12) and (13) were numerically solved via a finite difference method on a 2D 200 × 200 regular square lattice. Herein, space and time were discretised via central finite difference and forward Euler schemes, respectively. In order to solve eqn (12) and (13), both periodic and Neumann boundary conditions were employed, with the exact implementation varying between simulations.

For the validation of the PF model in Section 4.1, periodic boundary conditions are applied at all edges of the simulation domain. In the remaining sections, dealing with multiple aspects of the STrIPS process, periodic boundary conditions were only applied at the left and right (x̂) edges of the simulation domain. To control the diffusion of solvent, Neumann boundary conditions were used at the top and bottom (ŷ) edges instead. More detailed descriptions of these Neumann boundary conditions, along with their interpretations for the investigated systems, are provided in the respective Results sections and Section S0 of the ESI.†

As eqn (13) requires compositional inhomogeneities to induce phase separation, all simulations were started by seeding the initial composition ϕ₀ with random thermal noise, sampled uniformly from the [−0.01, 0.01] interval. Over the course of the resulting phase separation, the local interfacial tension σ was continuously calculated via eqn (8). When σ exceeded the critical value σ_c, phase separation was arrested by multiplying the value of the local relative mobility M̂ with a factor 10⁻⁶. The simulations were terminated when σ ≥ σ_c for the entirety of the simulation domain.

These respective values of σ_c were based on profiles of the average interfacial tension σ_av, showing the evolution of σ_av over the course of the simulations. The profiles of σ_av, along with the corresponding values of σ_c and a list of the other simulation parameters used in this work, can be found in the ESI† (Section S1).

Finally, it is worth to specifically mention that for all simulations the lattice spacing was kept at Δx̂ = Δŷ = 1 l̂. For a total system length of L = 100 μm, a typical length scale for the STrIPS system, this entails a characteristic length of l = 0.5 μm. With D_s = 10³ μm² s⁻¹ as a representative diffusion coefficient for a STrIPS solvent, this then gives a characteristic time of τ = 0.25 × 10⁻³ s.

4 Results

4.1 Model validation: construction of the phase diagram

First, the phase diagram in Fig. 2 was constructed to validate the PF model with respect to the underlying theory. In particular, the spinodal (dashed) and binodal (solid) curves were calculated in accordance with eqn (4) and (6) for χ₀ = 3.00 and ϕ_sc = 0.50 to qualitatively match the phase diagram of a typical STrIPS precursor mixture.^18,38 These theoretical curves delineate different regions of the phase diagram with distinct phase separation dynamics. PF simulations were then conducted for initial compositions across these different regions and the subsequent phase separation behaviour compared to theoretical expectations.

Fig. 2 The left phase diagram shows the spinodal (dashed) and binodal (solid) curves for STrIPS parameters χ₀ = 3.0 and ϕ_sc = 0.50. The markers give the initial compositions of performed simulations, while their shapes indicate the associated morphologies: spinodal channels (rhombus), elongated droplets (hexagon), dispersed circular droplets (pentagon), polydisperse nuclei (circle), homogeneous (square). The regions qualitatively corresponding to these morphologies are indicated in the right phase diagram.

For initial compositions under the spinodal curve, the results of the PF model are in full agreement with theoretical expectations. Namely, in this region of the phase diagram the oil/water mixture should be thermodynamically unstable and undergo spontaneous phase separation due to minor thermal fluctuations in the composition. The PF model readily captures this behaviour, known as spinodal decomposition, by displaying spontaneous demixing for initial compositions seeded with some uniformly sampled random noise (see Section 3.2).

Considering that this spinodal decomposition effectively provides the template for the STrIPS bijel structure, the influence of the initial composition on the PF system morphology is further investigated. Around the symmetry axis of the phase diagram (ϕ = 0.50), the PF simulations yield the characteristic spinodal morphology of continuous, interwoven channels (rhombus markers). However, moving away from the symmetry axis these continuous channels start to break up, forming elongated droplets of a minority phase in a continuous majority phase (pentagon markers). Consequently, in order for the STrIPS process to yield bicontinuous structures it is imperative that phase separation occurs within the centre of the spinodal region.

In direct contrast, for initial compositions close to the spinodal curve the PF model yields dispersions of distinctively circular droplets of the minority phase in a continuous majority phase (hexagon markers). It must be stressed that these droplets are formed due to spinodal decomposition, rather than the more commonly associated process of nucleation and growth. This notion is further supported through comparison with the PF results for initial compositions located in between the binodal and spinodal curve, where nucleation and growth is the dominant phase separation mechanism.

Namely, for initial compositions between the spinodal and binodal curves, the PF model fails to phase separate and returns homogeneous systems. This is in line with theory, as these compositions lie within the metastable region of the phase diagram. Herein, the oil/water mixture is stable with respect to small fluctuations in the composition and spinodal decomposition does not occur. However, phase separation can still happen through the nucleation of regions with a sufficiently different composition, which subsequently grow and coalesce. Such nuclei are not naturally formed by the PF model, which includes only minor random compositional noise, but can be introduced manually.

By seeding an initial composition in the metastable region with circular nuclei of the opposite (with respect to the symmetry axis) binodal composition, the PF simulations show the growth of these droplets in a manner analogous to Ostwald ripening, with larger droplets growing at the expense of the smaller ones. The resulting morphology resembles polydisperse circular droplets suspended in a continuous phase (circular markers). Consequently, while the PF model lacks the ‘nucleation’ component required for phase separation via the nucleation and growth mechanism, it adheres to theory by capturing the 'growth' aspect of such nuclei in the metastable region of the phase diagram.

Finally, for initial compositions outside the binodal curve PF simulations consistently return homogeneous systems, even after the manual addition of compositional nuclei. This is in agreement with theory, as for compositions outside the binodal curve the oil/water mixture should be unconditionally stable.

To summarise, the results of the PF model are fully consistent with the presented theory, showcasing different mechanisms of phase separation in their expected regions of the phase diagram. Even though the PF model can only partially reproduce the process of nucleation and growth, it fully describes phase separation via spinodal decomposition. Considering that the structure of the bijel mainly originates from spinodal decomposition, the PF model is therefore suitable for simulating the formation of bijels via STrIPS, as is done in the following.

4.2 Simulation of STrIPS bijel formation

With the theoretical soundness of the model firmly established, the next step was to examine whether the simplified, core components distilled from the STrIPS process are sufficient to emulate bijel formation. Specifically, it was investigated whether the PF model, operating under these basic guidelines, could produce structures with characteristic morphological features of the STrIPS bijel. To this end, the PF simulations needed to replicate the experimental conditions of a suitable STrIPS synthesis method.

Owing to its simple geometry and relative relevance for applications, the roll-to-roll fabrication of bijel films¹⁹ was chosen as the initial STrIPS model system. A schematic depiction of bijel formation in this system is shown in Fig. 3A. Here, the bijel precursor is deposited as a thin film on top of a solid substrate, which is immersed in a liquid ambient phase. The solid substrate blocks the diffusion of solvent from the bottom of the precursor film, causing phase separation to start only from the top interface with the ambient phase. As a result, bijel films formed in this manner display a characteristic trend in the average pore size: small pores at the top surface that progressively get larger closer to the solid substrate. This size gradient serves as an effective hallmark for the accuracy of the PF simulations, requiring the implementation of asymmetric diffusion to properly reflect the experimental conditions of the model system.

Fig. 3 (A) Schematic depiction of bijel film formation during roll-to-roll STrIPS. A precursor film is deposited on top of a solid substrate, while submerged in a liquid ambient phase. The solvent diffuses from the precursor mixture into the ambient phase, triggering phase separation and subsequent bijel formation. (B) PF simulations of bijel film formation in the roll-to-roll model system. Starting from a precursor mixture with ϕ₀ = 0.50 ± 0.01 and ϕ_s = 0.50, the critical composition for values of χ₀ = 3.0 and ϕ_sc = 0.50, the diffusion of solvent into the ambient phase causes the interaction parameter χ to increase (bottom row). This induces phase separation in the oil/water mixture ϕ (top row), in turn giving rise to an increasing interfacial tension σ (central row). Where this σ exceeds the critical values σ_c = 1.0, phase separation is locally arrested. After t̂ = 1500 (t = 3.75 s for a film thickness L = 100 μm), phase separation and structure coarsening has ceased over the entire simulation domain and the simulation is terminated.

This asymmetric diffusion of solvent is readily incorporated into the PF model through different boundary conditions. Based on the schematic depiction of the model system in Fig. 3A, the simulation domain represents a 2D cross-section of the deposited precursor film, made to be infinitely wide via periodic boundary conditions at the left and right edges. Employing a Neumann boundary condition, solvent is allowed to diffuse across the top edge into the ambient liquid. Considering that the volume of the ambient phase is significantly larger than that of the precursor film, the solvent content of the ambient phase ϕ^BC_s is treated as constant. At the same time, another Neumann boundary condition prevents solvent from crossing the bottom edge of the system, representing the presence of a solid substrate.

Consistent with the different interpretations of ϕ and ϕ_s in Section 2.2, exchange between the precursor and the ambient liquid is thus limited solely to the solvent fraction ϕ_s. This is further enforced by imposing Neumann boundary conditions on eqn (13) that close both the top and bottom system boundary to the direct transport of oil/water associated with phase separation. With these different boundary conditions, PF simulations were conducted that qualitatively matched the experimental conditions of the STrIPS roll-to-roll process, the results of which are depicted in Fig. 3B.

Specifically, the images in Fig. 3B illustrate the different stages of structure formation in the PF simulations of STrIPS. Phase separation is shown through the evolution of the oil fraction ϕ, along with the associated local interfacial tension σ. The diffusion of the solvent is represented by the oil/water interaction parameter χ. Mirroring the experimental procedure, the simulation starts (t̂ = 0) from a homogeneous precursor with a composition located at the critical point of the phase diagram. For the aforementioned STrIPS precursor parameters of χ₀ = 3.0 and ϕ_sc = 0.50, this corresponds to an initial composition of ϕ₀ = 0.50 ± 0.01 and ϕ_s = 0.50 (see Fig. 1). As a homogeneous mixture, the initial state has constant values for χ and σ, with notably σ = 0. This changes, however, as the difference in solvent content between the precursor mixture and the liquid ambient phase, which contains no solvent (ϕ^BC_s = 0), causes the solvent to diffuse out of the precursor into the ambient phase.

Namely, the diffusion of solvent shifts the composition in the adjacent region of the precursor mixture into the spinodal region. The corresponding increase in χ induces local phase separation (t̂ = 3750), with the resulting formation of interfaces giving rise to a non-zero value for σ. As the solvent continues to diffuse away, this front of phase separation and increasing interfacial tension progressively makes its way through the system (t̂ = 7500). At a certain point, σ locally exceeds the set threshold value, here σ_c = 1.0, resulting in the arrest of phase separation in that region. Meanwhile, phase separation and structure coarsening continue in the deeper regions where the interfacial tension remains lower. Eventually, the solvent level will have been sufficiently depleted to induce phase separation over the entire system domain (t̂ = 11 250), with the evolution of the interfacial tension σ lagging somewhat behind. Finally, σ catches up and causes the complete arrest of phase separation and structure coarsening throughout the entire simulation domain (t̂ = 15 000), after which the simulation is stopped.

These simulations evidently demonstrate that the PF model can capture the key components of the STrIPS process. Not only do they display progressive phase separation induced by the diffusion of solvent, they also indicate that the interfacial tension σ functions as an effective measure for arresting phase separation. The most obvious sign, however, is the fact that the resulting morphology bears a striking resemblance to the STrIPS bijel as expected for the roll-to-roll process.

Most notably, the final morphology in the top row of Fig. 3B has a distinct gradient in the pore size. Matching experimental results, the pores are smallest at the top edge of the system domain, at the interface with the ambient liquid. Deeper into the structure the pores get gradually larger, reaching a maximum at the bottom next to the solid substrate. As this characteristic morphological feature of the STrIPS bijel was reproduced by the PF model, its origin in simulations could help elucidate their appearance in experiments as well.

In the PF simulations, the pore size gradient results due to there being local differences in the solvent flux over the course of the simulation. Deeper in the structure, the solvent content remains higher, for longer, compared to the region adjacent to the ambient phase. Consequently, there is a longer period of time between the onset of phase separation and the interfacial tension reaching σ_c, arresting the structure. As the structure continues to coarsen in this period, this results in comparatively larger pores than near the surface. Notably, these observations support the proposed explanations of pore size gradients for experimental STrIPS bijels.^8,19

While these findings highlight the potential of the PF model for simulating bijel formation via STrIPS, some major points require further discussion. In particular, the influence of the PF model parameters on the obtained morphology and how this physically relates to the STrIPS process. For example, the simulations in Fig. 3 were performed for M̂ = 0.01 and [small kappa, Greek, circumflex] = 0.50 in eqn (13). The values of these parameters and their respective interpretations proved to be crucial for the success of the PF model in simulating STrIPS. As such, they are discussed in detail in the next section.

4.3 The influence of M̂ and [small kappa, Greek, circumflex]

The dynamics of phase separation and structure coarsening are crucial elements of the formation of bijels via STrIPS. In this PF model, these elements are represented through the relative oil/water mobility M̂ and gradient energy coefficient [small kappa, Greek, circumflex] in eqn (13). To gain further insight into the influence of these parameters on the formed morphology, STrIPS simulations as in Fig. 3 were performed with different combinations of M̂ and [small kappa, Greek, circumflex]. In particular, somewhat low values of M̂ and [small kappa, Greek, circumflex] were expected to yield the characteristic spinodal morphology associated with the bijel.

Namely, these low parameter values represent a system where solvent diffusion is fast compared to phase separation (M̂) and the presence of compositional gradients is not severely penalised ([small kappa, Greek, circumflex]). The former effectively imposes a ‘deep quench’ on the bijel precursor, bringing the composition deeper into the spinodal region prior to phase separation, whereas the latter reduces the influence of local composition inhomogeneities. Reflecting the STrIPS process, this combination should facilitate spinodal decomposition in an isotropic manner over the entire depth of the system.

This expectation is seemingly confirmed by the results of the PF model, shown in Fig. 4 for different combinations of M̂ and [small kappa, Greek, circumflex]. Lower values of M̂ and [small kappa, Greek, circumflex] result in the formation of structures visually similar to experimental STrIPS bijels, containing randomly oriented domains with a size gradient perpendicular to the surface. In contrast, higher values of M̂ and [small kappa, Greek, circumflex] seem to favour the formation of a more laminar morphology after phase separation, displaying alternating layers of both the oil- and water-rich phase parallel to the surface. Close to the interface with the ambient phase, however, these laminar structures contain a localised region with a notably different structure. This surface region displays a wider variety of domain orientations, appearing similar to the structures observed for lower M̂ and [small kappa, Greek, circumflex]. Overall, however, it can be concluded that low values of M̂ and [small kappa, Greek, circumflex] yield isotropic bijel structures, whereas higher values induce the formation of laminar pores.

Fig. 4 Final morphologies of the phase-separated oil/water mixtures after STrIPS for different values of the relative mobility M̂ and the gradient energy coefficient [small kappa, Greek, circumflex]. Lower values of M̂ and [small kappa, Greek, circumflex] yield relatively isotropic phase separation, whereas higher values result in laminar structures with a more isotropic surface region.

As illustrated by simulations shown in Fig. 5, this observed morphological differences can be explained by progressive phase separation directed by a planar diffusion front. That is, phase separation sequentially occurring in distinct, adjacent regions of the simulation domain, following the sufficient depletion of solvent. In particular, it is the minimisation of the gradient energy associated with this mechanism of phase separation that causes the morphological differences. Namely, compositional differences between adjacent regions will induce the formation of a boundary layer in order to minimise the gradient energy. For a planar diffusion front, these boundary layers can give rise to the observed laminar structures. Whether this actually occurs, however, depends on the extent of progressive phase separation.

Fig. 5 Progressive phase separation during STrIPS for low (i) and high (ii) values of the relative mobility M̂ and gradient energy coefficient [small kappa, Greek, circumflex]. For low values (M̂ = 0.01 and [small kappa, Greek, circumflex] = 0.50), phase separation can not keep up with the solvent diffusion front, qualitatively indicated by the white dashed line. This results in the relatively isotropic phase separation typical for STrIPS bijels. For high values (M̂ = 0.05 and [small kappa, Greek, circumflex] = 1.00), phase separation can keep up with the diffusion front, causing the formation of planar boundary layers that subsequently coarsen into laminar pores.

As it controls the relative dynamics of phase separation and solvent diffusion, this extent is mainly set by M̂. High values of M̂ reflect systems with relatively fast phase separation compared to solvent diffusion. In these systems, phase separation can keep up with the moving diffusion front, causing the edge of this diffusion front to effectively become an interface between the region undergoing phase separation and the adjacent region where it has yet to occur. Minimisation of gradient energies will then induce the formation of a planar boundary layer, matching the shape of the diffusion front (ESI,† S3). As the diffusion front progressively moves deeper into the structure, the associated phase separation continuously generates such planar boundary layers, which eventually coarsen into the laminar morphologies observed in Fig. 4 and 5.

In contrast, for low values of M̂ the dynamics of phase separation are slow compared to the solvent diffusion. In such systems, phase separation can not keep up with the moving diffusion front, thereby lowering the influence of the associated interface and allowing for phase separation in a more isotropic manner. The same principle explains the existence of the differently structured surface regions observed for higher values of M̂. As these regions are directly adjacent to the ambient phase, phase separation can not keep up with the rapid depletion of solvent from this region. This rapid depletion of solvent causes the formation of a locally isotropic structure, which is then quickly arrested before considerable coarsening can take place.

The influence of [small kappa, Greek, circumflex] on the morphology is similar to M̂, albeit less pronounced. The value of [small kappa, Greek, circumflex] determines the relative energy penalty for gradients in the composition, thereby also setting the degree of influence a moving interface such as the diffusion front can exert over the system. For higher values of [small kappa, Greek, circumflex] smaller compositional variations become more relevant, enhancing the ability of the diffusion front to direct phase separation. This is mainly noticeable at lower values of M̂, such as in the left column of Fig. 4. Herein, the increase in [small kappa, Greek, circumflex] causes the formation of more laminar structures, mostly located in the central region of the structure. Both the regions located near the surface (top) and the solid support (bottom), however, display relatively isotropic features. This is because progressive phase separation is most prominent in the central region. As mentioned above, the surface region has a more isotropic structure due to the rapid depletion of solvent. For the region directly adjacent to the solid surface, the shape of solvent profile remains mostly flat over time (ESI,† S3). As such, phase separation occurs relatively simultaneously throughout this region and is therefore not significantly influenced by a moving diffusion front, also yielding a more isotropic morphology.

Concluding, the PF simulations revealed that the interplay between phase separation and structure coarsening, represented by the values of M̂ and [small kappa, Greek, circumflex], plays a major role in acquiring isotropic bijels via STrIPS. In particular, it was found that isotropic bijel structures only occurred for fast solvent diffusion compared to phase separation, while different relative dynamics resulted in the formation of laminar pores.

With these relations established, the final parameter in the PF model to be discussed is the critical interfacial tension σ_c. Similar to M̂ and [small kappa, Greek, circumflex], the value of σ_c has a distinct influence on the morphology and can be directly related to a physical property of the bijel system, as will be shown next.

4.4 The influence of σ_c

Owing to its considerable influence on the obtained morphology, the value of the critical interfacial tension σ_c warrants additional discussion. The relevance of σ_c, above which phase separation is arrested, can be readily made apparent by considering it as an effective stand-in for the nanoparticle size. Namely, in accordance with eqn (14), the attachment energy ΔG of a particle at a liquid–liquid interface scales with the interfacial tension σ times the square of the particle radius r²³⁹

ΔG = −πr²σ(1 − |cos θ|)² (14)

with additionally θ as the contact angle of the nanoparticle at the liquid interface.

Consequently, as an interface develops during phase separation, the increasing interfacial tension will irreversibly trap larger particles earlier than smaller ones. In other words, larger particles can cause interfacial jamming at lower interfacial tensions compared to smaller particles. This can be qualitatively reflected in the value of σ_c, with lower values corresponding to the presence of larger particles in the system. In this manner, the pore size can be controlled via σ_c, with lower values of σ_c resulting in smaller pores.

To investigate this potential effect of σ_c, PF simulations were carried out with different values of σ_c, while keeping the other parameters the same as for the STrIPS simulations in Fig. 3. Some representative examples of the resulting morphologies are shown in Fig. 6A, all displaying the characteristic STrIPS bijel structure with a distinct pore size gradient. While the pore sizes visually already seem to differ between them, the pore size profiles in Fig. 6B were calculated (ESI,† S4) to properly compare the different morphologies.

Fig. 6 (A) Representative morphologies obtained after STrIPS PF simulations for different values of the critical interfacial tension σ_c. (B) Profiles of the average pore size d̂/L̂ over the depth of the simulation domain ŷ/L̂ for different values of the critical interfacial tension σ_c. Note that both the pore size and the depth have been normalised with respect to the simulation domain length L̂.

These profiles show the average pore size over the depth of the structure, both normalised with respect to the length of the simulation domain, for different values of σ_c. Corroborating previous results, they unambiguously demonstrate that the average pore size increases deeper into the structure, further away from the boundary with the ambient phase. Furthermore, the profiles also indicate that the value of σ_c does have a noticeable effect on the pore size. In line with expectations, lower values of σ_c translate to smaller average pore sizes over the entire depth of the structure. The strength of this effect, however, varies within the structure, being the most pronounced close to the surface and lessening deeper in.

This difference can be explained by considering that the surface region, being directly adjacent to the ambient phase, is depleted of solvent considerably faster compared to the rest of the system. Therefore, the quench of the composition into the spinodal region is faster and deeper, resulting in correspondingly fast dynamics of phase separation and a relatively steep increase in the interfacial tension. As such, equivalent changes in σ_c result in greater differences in the average pore size for the structure's surface compared to ‘slower’ regions deeper within. Regardless of this discrepancy, it can be concluded that σ_c can effectively be used to tune the particle size in the bijel system.

While these findings seem to align with the expected behaviour of σ_c in its qualitative relation to the particle size, there remain some key issues with directly transposing it to the experimental system. Most importantly, the PF model does not sufficiently take into account the explicit diffusion, attachment and subsequent jamming of particles.

For example, not only do smaller particles simply diffuse faster, the bijel structure is also not immediately arrested following the irreversible attachment of particles, as implied by the σ_c rationale in the PF simulation. Rather, structure coarsening continues until the decreasing interfacial area causes the jamming of an attached particle layer. Consequently, it is the liquid interfacial area available, in relation with the number of particles present, that sets the criterion for arresting phase separation. For example, when present in equivalent volumes, smaller particles require larger interfacial areas compared to larger particles, resulting in the stabilisation of smaller pores.^15,40 Furthermore, the minimum dimensions of a particle-stabilised pore are dictated by the size of the particles. That is, particles can not stabilise pores with smaller dimensions than the particles themselves. Consequently, for large differences in particle size, the smaller particles are likely to yield structures with smaller pores compared to larger ones.

Even though its relation with the physical bijel might not be fully complete, the value of σ_c remains an effective tool for controlling the average pore size in PF simulations. With the role of this parameter elucidated, in addition to M̂ and [small kappa, Greek, circumflex], the next step was to directly compare PF simulations to experimental bijel systems. In particular, it is demonstrated that the PF model can match trends in bijel morphologies for different experimental conditions, providing a potential explanation behind their origin.

4.5 Comparison with experimental results

So far, the PF model has shown capable of replicating the core components of the STrIPS process, yielding structures with very similar morphological features to bijels. To further support this notion, in addition to demonstrating the potential of the PF model in rationalising experimental findings, PF simulations are directly compared to experimental morphologies of STrIPS bijels synthesised under different experimental conditions. As will be shown, the observed morphological trends can be readily explained by the PF model, providing further insight into the tailoring of the bijel structure for this experimental system.

The experimental system in questions involves the formation of a bijel fibre containing a central borehole, which is filled with a liquid. While the experimental details of this system are described elsewhere,⁴¹ a schematic depiction of the PF model is shown in Fig. 7A. The principal feature of this specific bijel system is that the solvent can diffuse in two directions during STrIPS, both into the exterior ambient phase and into the liquid occupying the central bore. In experiments, it is found that this results in the formation of a two-sided pore size gradient: the pores are smallest next to the interfaces with both the interior bore liquid and the ambient phase, while they are largest near the centre of the bijel shell.

Fig. 7 (A) Schematic depiction of the fabrication of bijel fibres containing a bore channel via STrIPS. A precursor mixture and a bore liquid are coaxially extruded into a liquid ambient phase, resulting in the formation of cylindrical bijel shells. Shown in the inset, the PF model represents a partial cross-section of this cylindrical shell as a square domain, neglecting curvature. Here, solvent is allowed to diffuse both into the bore liquid (bottom) and the ambient phase (top). (B) STrIPS PF simulations of the model system shown in (A), using different values of the solvent fraction in the bore liquid ϕ^BC_s. For increasing ϕ^BC_s, solvent diffusion and phase separation become increasingly asymmetric, influencing the pore size gradient present in the final morphology.

Additionally, the solvent content in the bore liquid can be varied. Experimentally, it was found that varying this solvent content shifts the position of the maximum pore size in the bijel shell, moving closer to the interface with the bore liquid for higher solvent contents. In order for simulations to assist in providing a possible explanation behind the origin of these morphological trends, the experimental system first had to be replicated in the PF model.

Fig. 7A shows a schematic depiction of the PF model used for this experimental system. A section of the precursor shell is represented by a square simulation domain, neglecting the influence of curvature on the diffusion of solvent. The top edge of the simulation domain represents the boundary between the bijel precursor and the exterior ambient phase, whereas the bottom edge represents the boundary with the interior bore channel. The exchange of solvent with both the liquid bore channel and the ambient phase was controlled through Neumann boundary conditions. In particular, the solvent content in the liquid bore channel was varied through ϕ^BC_s, whilst for the ambient phase it was consistently kept at zero. Subsequently, PF simulations were performed with the same parameters as for the STrIPS simulations in Fig. 3, while the values of ϕ^BC_s were varied.

The results of these simulations, shown in Fig. 7B, reveal that the diffusion of solvent into the ambient phase and liquid bore channel initiates phase separation at the top and bottom of the simulation domain, respectively. When there is no solvent present in the bore liquid (ϕ^BC_s = 0), this phase separation continues in a symmetric fashion towards the centre of the system. As the solvent content in the centre remains higher over a longer period of time, prolonged coarsening of the structure results in larger pores. This process finally yields a morphology with a distinct two-sided pore size gradient, corroborating experimental observations.

The introduction of solvent into the bore liquid, however, induces asymmetric phase separation that considerably changes the resulting morphology. In particular, higher solvent contents in the bore liquid lower the rate of solvent depletion in the adjacent precursor region, delaying the onset of phase separation. Additionally, this region maintains higher levels of solvent over a longer time, allowing for a greater extent of coarsening prior to the arrest of phase separation. Consequently, the pores closer to the bore channel become larger with increasing ϕ^BC_s. An overview of the resulting structures for the different values of ϕ^BC_s is shown in Fig. 8A, additionally containing the corresponding morphologies of experimental STrIPS bijels, acquired via confocal microscopy.

Fig. 8 (A) Partial cross-sections of bore-channel bijels for different solvent fractions in the interior bore liquid ϕ^BC_s, obtained both experimentally (confocal microscopy) and via PF simulations. The structure depth within the experimental and simulated cross-sections is normalised with respect to length L̂, corresponding to the length of the confocal image and simulation domain, respectively. In this orientation, the interfaces with the ambient phase and bore channel are respectively located at ŷ/L̂ = 0 and ŷ/L̂ = 1. (B) Pore size profiles over the experimental and simulated morphologies of the bore-channel bijel, obtained for different values of the bore liquid solvent fraction ϕ^BC_s. The average pore size d̂/L̂ is calculated over the depth of the structure ŷ/L̂ and the resulting points fitted with a quadratic polynomial. The position of the maximum pore size is indicated by the star marker, shifting with increasing ϕ^BC_s. Note that the legend provided for the simulated morphologies (left) also applies to the experimental morphologies (right).

On first sight, the experimental and simulated morphologies of the STrIPS bijels already appear strikingly similar, indicating that the PF model accurately managed to replicate the experimental conditions of the STrIPS synthesis. To quantitatively substantiate this observation, pore size profiles were calculated for the different values of ϕ^BC_s and subsequently fitted with a quadratic polynomial curve. In order to compare the pore size profiles between the simulations and the experimental images, the average pore size d̂ and the depth within the structure ŷ were normalised with respect to the length L̂, corresponding to the length of the simulation domain and confocal image, respectively. These profiles, shown in Fig. 8B, confirm the presence of a two-sided pore size gradient for both the PF simulations and the experimental bijel morphologies. Starting from either the top or the bottom boundary of the structure the pores are the smallest, which then increase in size deeper into the structure before reaching a maximum and subsequently decreasing again while approaching the opposite boundary. The depth where the pore size reaches its maximum, indicated by the star marker in the pore size profiles, appears to shift depending on the solvent fraction of the bore liquid.

To better visualise this shift for both experimental and simulated morphologies, Fig. 9 shows the relative depth of maximum pore size ŷ_max/L̂ plotted against the solvent volume fraction within the bore liquid ϕ^BC_s. As seen in the simulations of Fig. 7B, phase separation for ϕ^BC_s = 0 is relatively symmetric due to the equal rates of solvent depletion to the ambient phase and bore channel. This results in an equally symmetric, parabolic pore size profile with a maximum located near the centre of the morphology, roughly around ŷ/L̂ = 0.5. In contrast, increasing the solvent fraction in the bore liquid slows down the diffusion rate towards the bore channel relative to the ambient phase. As a result, phase separation in the region adjacent to the bore channel is arrested at a later stage, facilitating a greater extent of structure coarsening. Consequently, the average pore size on the side of the bore channel increases, shifting the maximum away from the centre towards the interface with the bore liquid, corresponding to an increase in ŷ_max/L̂.

Fig. 9 The maximum pore depth ŷ_max/L̂, corresponding to the relative depth within the bijel cross-section where the pore size reaches its maximum value, versus the solvent fraction of the interior bore channel ϕ^BC_s. The green markers correspond to experimental bijel structures, whereas the blue markers correspond to the structures provided by the phase-field model.

While this shift in the location of the maximum pore size as a function of ϕ^BC_s is visible for both the experimental and simulated morphologies, it is distinctly more pronounced in the simulations. There, increasing the solvent content in the bore liquid up to ϕ^BC_s = 0.40 nearly creates a linear pore size profile with a maximum very close to the boundary of the bore channel. For the experimental structures, however, the shift in the pore size maximum appears to flatten off somewhat, not exceeding ŷ_max/L̂ = 0.63 for ϕ^BC_s = 0.40.

Likely, part of the discrepancy can be attributed to the fact that the solvent parameter ϕ_s in simulations simply does not directly translate to the actual solvent volume fraction as employed in experimental work. Additionally, there are other, underlying mechanisms at play during STrIPS that are currently not captured by the model. For example, the influx of both ambient and bore liquids or the preferential partitioning of solvents¹⁸ can disrupt the ideal diffusion behaviour as assumed by the model, causing deviation from the predicted pore size profiles.

In conclusion, there is a significant level of agreement between the simulated and experimental bijel morphologies, indicating that the PF model can accurately describe bijel formation via STrIPS in this complex experimental system. In this manner, the PF simulations allowed for the in situ observation of the STrIPS process, providing a potential explanation for the morphological trends found in the experimental bijels depending on the synthesis conditions.

5 Discussion

The results presented in this work demonstrate that the simple PF model can effectively capture the core mechanisms of the complex STrIPS process. Through these basic principles, the simulations readily produce structures with the characteristic features of STrIPS bijels. Consequently, the PF model allows for the observation and investigation of the entirety of the STrIPS process, not only the final structure. In this manner, the influence of critical parameters, such as the mobility M̂, on the final morphology could be directly associated with differences in the formation mechanism. Moreover, the PF simulations were used to rationalise the appearance of morphological trends in experimental STrIPS bijels.

While these findings showcase the potential of the PF model, its current version remains limited in some aspects. Firstly, the approach of modelling phase separation and solvent diffusion through two separate phase-fields prevents direct comparison between the model solvent parameter ϕ_s and the experimental solvent volume fraction. Additionally, this method excludes multi-component interactions and by extension the associated effects that could be relevant for the STrIPS process, such as preferential solvent partitioning.¹⁸

Secondly, the PF model employs considerable simplifications by assuming constant values of the mobility M and a free energy density f given by a regular solution. While these are sufficient for describing the main principles underlying STrIPS, they need to be further extended to provide a more accurate representation. For example, composition-dependent mobilities can result in different coarsening kinetics for the formed structure.⁴² Additionally, a Flory–Huggins free energy density could better match the thermodynamics of the STrIPS precursor mixture while simultaneously including multi-component interactions.³⁸ This extension of the current model to include full multi-component phase separation with composition-dependent mobilities will be the topic of future work.

Furthermore, the present version of the PF model only includes diffusive mass transport, neglecting convective transport via fluid flows. However, as is the case for polymeric membranes fabricated via phase inversion,^30,33 it is likely that such hydrodynamic effects also play a role in determining the final morphology of the STrIPS bijel. Consequently, a full description of the STrIPS process requires the incorporation of hydrodynamics. While outside the scope of this current work, the coupling of hydrodynamics into the PF model has the additional benefit of expanding the selection of relevant system parameters that can be examined. For example, the influence of the component viscosity on the morphology of the STrIPS bijel, currently reflected only to a limited extent in the relative mobility M̂, could be explicitly investigated.

Finally, while the interfacial jamming of nanoparticles during STrIPS is effectively replicated through the surface tension σ, there remain some issues with this approach (see Section 4.4). PF models, however, have additional options for the direct inclusion of the nanoparticles. For example, ongoing research involves describing the nanoparticles through their own continuous phase-field. In this manner, the interfacial activity and subsequent arrest of phase separation could be modelled while remaining computationally feasible on the spatial and temporal scales relevant to STrIPS.

6 Conclusions

In this work, a PF model is presented that captures the complex mechanism of bijel formation via STrIPS through relatively simple methods. This is achieved through the evolution of two distinct phase-fields, representing the oil/water mixture and the solvent, respectively. By making the oil/water interaction dependent on the solvent field, the physics of the oil/water/solvent system can be effectively replicated, showing phase separation behaviour consistent with the underlying theory. It is further demonstrated that the role of the nanoparticles in arresting phase separation can be reproduced by using the interfacial tension as a substitute measure.

Through this PF model, it is then shown that morphologies with characteristic bijel features are obtained in systems with relatively slow dynamics of phase separation compared to solvent diffusion. In the model, such a system is represented by low values of the relative oil/water mobility M̂ and gradient energy coefficient [small kappa, Greek, circumflex]. Increasing these values, thus transitioning to a system with faster relative dynamics and a higher energy penalty for compositional gradients, instead results in the formation of laminar pores. These morphological differences are explained in the context of progressive phase separation directed by a planar diffusion front.

Finally, the PF model is applied to investigate the formation of STrIPS bijels by replicating the conditions of a complex experimental system. The resulting, simulated morphologies show strong qualitative agreement with the experimental bijel structures obtained from experiments. Moreover, the PF model could reproduce morphological trends found in the experimental system, while simultaneously providing a potential explanation for their origin. While some discrepancies remain, it is expected that these can be further remedied by extending the model to include full multi-component phase separation.

Overall, this work unambiguously demonstrates the potential of the PF model in the simulation of bijel formation via STrIPS, especially in conjunction with experimental work. By allowing effective in situ observation of the entire STrIPS process, the PF model can help clarify the formation mechanisms of current STrIPS bijel morphologies, as well as facilitate the synthesis of new ones.

Author contributions

Jesse Steenhoff: conceptualisation, methodology, software, investigation, formal analysis, data curation, visualisation, writing – original draft, writing – review & editing Martin Haase: conceptualisation, supervision, funding acquisition, writing – review & editing.

Data availability

The code and datasets underlying the findings in this work are available on GitHub (https://github.com/JesseSteenhoff/A-simple-phase-field-model-of-bijel-formation-via-solvent-transfer-induced-phase-separation) and have been archived on Zenodo (https://doi.org/10.5281/zenodo.14755856). These include the Python scripts for the presented phase-field model and the pore size analysis. The datasets (both experimental and simulated) required for the reproduction of the pore size analysis are additionally provided.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This publication is part of the project “Bijel templated membranes for molecular separations” (with project number 18632 of the research programme Vidi 2019), which is financed by the Dutch Research Council (NWO). The authors would like to acknowledge Mariska de Ruiter for providing the raw confocal microscopy data of the bijel structures with bore channels.

References

1. R. Aveyard, B. P. Binks and J. H. Clint, Adv. Colloid Interface Sci., 2003, 100–102, 503–546.
2. Y. Yang, Z. Fang, X. Chen, W. Zhang, Y. Xie, Y. Chen, Z. Liu and W. Yuan, Front. Pharmacol., 2017, 8, 235054.
3. H. Zhao, Y. Yang, Y. Chen, J. Li, L. Wang and C. Li, Chem. Eng. Sci., 2022, 248, 117085.
4. J. Texter, Colloid Polym. Sci., 2022, 300, 587–592.
5. K. Stratford, R. Adhikari, I. Pagonabarraga, J.-C. Desplat and M. E. Cates, Science, 2005, 309, 2198–2201.
6. E. M. Herzig, K. A. White, A. B. Schofield, W. C. K. Poon and P. S. Clegg, Nat. Mater., 2007, 6, 966–971.
7. M. F. Haase, H. Jeon, N. Hough, J. H. Kim, K. J. Stebe and D. Lee, Nat. Commun., 2017, 8, 1234.
8. H. Siegel, A. J. Sprockel, M. S. Schwenger, J. M. Steenhoff, I. Achterhuis, W. M. de Vos and M. F. Haase, ACS Appl. Mater. Interfaces, 2022, 14, 43195–43206.
9. T. Wang, R. A. Riggleman, D. Lee and K. J. Stebe, Mater. Horiz., 2023, 10, 1385–1391.
10. J. A. Witt, D. R. Mumm and A. Mohraz, J. Mater. Chem. A, 2016, 4, 1000–1007.
11. D. Cai, F. H. Richter, J. H. J. Thijssen, P. G. Bruce and P. S. Clegg, Mater. Horiz., 2018, 5, 499–505.
12. G. Di Vitantonio, T. Wang, M. F. Haase, K. J. Stebe and D. Lee, ACS Nano, 2019, 13, 26–31.
13. T. Park, G. H. Choi, D. Lee and P. J. Yoo, J. Membr. Sci., 2021, 634, 119386.
14. J. W. Tavacoli, J. H. J. Thijssen, A. B. Schofield and P. S. Clegg, Adv. Funct. Mater., 2011, 21, 2020–2027.
15. D. Cai, P. S. Clegg, T. Li, K. A. Rumble and J. W. Tavacoli, Soft Matter, 2017, 13, 4824–4829.
16. M. F. Haase, K. J. Stebe and D. Lee, Adv. Mater., 2015, 27, 7065–7071.
17. S. P. Kharal and M. F. Haase, Small, 2022, 18, 2106826.
18. M. A. Khan, A. J. Sprockel, K. A. Macmillan, M. T. Alting, S. P. Kharal, S. Boakye-Ansah and M. F. Haase, Adv. Mater., 2022, 34, 2109547.
19. H. Siegel, M. de Ruiter, G. Athanasiou, C. M. Hesseling and M. F. Haase, Adv. Mater. Technol., 2024, 9, 2301525.
20. L. Bai, J. W. Fruehwirth, X. Cheng and C. W. Macosko, Soft Matter, 2015, 11, 5282–5293.
21. D. Cai and P. S. Clegg, Chem. Commun., 2015, 51, 16984–16987.
22. E. Kim, K. Stratford, R. Adhikari and M. E. Cates, Langmuir, 2008, 24, 6549–6556.
23. P. C. Millett and Y. U. Wang, J. Colloid Interface Sci., 2011, 353, 46–51.
24. F. Jansen and J. Harting, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2011, 83, 046707.
25. J. M. Carmack and P. C. Millett, J. Chem. Phys., 2015, 143, 154701.
26. E. Kim, K. Stratford and M. E. Cates, Langmuir, 2010, 26, 7928–7936.
27. A. S. Joshi and Y. Sun, Phys. Rev. E: Stat., Nonlinear, Soft Matter Phys., 2009, 79, 066703.
28. N. Karthikeyan and U. D. Schiller, Soft. Matter., 2024, 20, 8952–8967.
29. M. J. A. Hore and M. Laradji, J. Chem. Phys., 2007, 126, 244903.
30. B. Zhou and A. C. Powell, J. Membr. Sci., 2006, 268, 150–164.
31. Y. Mino, T. Ishigami, Y. Kagawa and H. Matsuyama, J. Membr. Sci., 2015, 483, 104–111.
32. G. Zhang, T. Yang, S. Yang and Y. Wang, Phys. Rev. E, 2017, 96, 032501.
33. D. R. Tree, T. Iwama, K. T. Delaney, J. Lee and G. H. Fredrickson, ACS Macro Lett., 2018, 7, 582–586.
34. M. R. Cervellere, X. Qian, D. M. Ford, C. Carbrello, S. Giglia and P. C. Millett, J. Membr. Sci., 2021, 619, 118779.
35. J. M. Carmack and P. C. Millett, Soft Matter, 2018, 14, 4344–4354.
36. J. W. Cahn and J. E. Hilliard, J. Chem. Phys., 1958, 28, 258–267.
37. R. C. Desai and R. Kapral, in Free energy functional, Cambridge University Press, Cambridge, 2009, pp. 20–24.
38. H. Zhang, Y. Wu, F. Wang, F. Guo and B. Nestler, Langmuir, 2021, 37, 5275–5281.
39. B. P. Binks and T. S. Horozov, Colloidal Particles at Liquid Interfaces, Cambridge University Press, Cambridge, 2006.
40. M. E. Cates and P. S. Clegg, Soft Matter, 2008, 4, 2132–2138.
41. R. de Ruiter, PhD thesis (Research UU/Graduation UU), Universiteit Utrecht, 2024.
42. B. König, O. J. J. Ronsin and J. Harting, Phys. Chem. Chem. Phys., 2021, 23, 24823–24833.

---

Footnote

† Electronic supplementary information (ESI) available. See DOI: https://doi.org/10.1039/d4cp04638b

---