How antagonistic salts cause nematic ordering and behave like diblock copolymers

TL;DR

This paper uses simulations and theory to show how antagonistic salts induce nematic ordering and complex structures in binary fluids, linking electrostatic interactions to the formation of lamellar phases and their dynamics.

Contribution

It introduces a theoretical framework connecting antagonistic salt effects to the Ohta-Kawasaki model, predicting structure sizes and morphologies in binary fluid mixtures.

Findings

Spinodal decomposition is arrested with small electrostatic interactions.

Lamellar structures exhibit nematic order and grow over time following a salt-dependent power law.

High salt concentrations lead to structure dissolution due to vanishing surface tension.

Abstract

We present simulation results and an explanatory theory on how antagonistic salts affect the spinodal decomposition of binary fluid mixtures. We find that spinodal decomposition is arrested and complex structures form only when electrostatic ion-ion interactions are small. In this case fluid and ion concentrations couple and the charge field can be approximated as a polynomial function of the relative fluid concentrations alone. When the solvation energy associated with transfering an ion from one fluid phase to the other is of the order of a few k_BT, the coupled fluid and charge fields evolve according to the Ohta-Kawasaki free energy functional. This allows us to accurately predict structure sizes and reduce the parameter space to two dimensionless numbers. The lamellar structures induced by the presence of antagonistic salt in our simulations exhibit a high degree of nematic ordering and the growth of ordered domains over time seems to follow a power law. This power law carries a time exponent proportional to the salt concentration. We reproduce and interpret neutron scattering data from previous experiments of similar systems. The dissolution of structures at high salt concentrations observed in these experiments agrees with our simulations and we explain it as the result of a vanishing surface tension due to electrostatic contributions. We conclude by presenting preliminary 3D results showing the same morphologies as predicted by the Ohta-Kawasaki model as a function of volume fraction and suggesting that our findings from 2D systems remain valid in 3D.

Equations

$$f_i⟶f_i-\frac{1}{\tau }(f_i-f_i^{\mathrm{eq}})+S_i.$$

$$\begin{array}{r}{\mathcal{F}}_T={\mathcal{F}}_F+{\mathcal{F}}_I=\stackrel{{\mathcal{F}}_F}{\overbrace{\int g(\psi ,\nabla \psi )}}+\\ \underset{{\mathcal{F}}_I}{\underbrace{\sum _{\pm }{k}_{B}T{c}_{\pm }[\ln (v_0{c}_{\pm })-1]+\sum _{\pm }{c}_{\pm }{\mu }_{\pm }^{\mathrm{sol}}(\psi )+\frac{1}{2}\rho \varphi \mathrm{d}V}}.\end{array}$$

$$\frac{\partial c_{\pm }}{\partial t}$$

$${\vec{J}}_{\pm }$$

$${\mu }_{c_{\pm }}=\frac{\delta {\mathcal{F}}_I}{\delta c_{\pm }}=k_{B}T\ln (v_0{c}_{\pm })+\frac{\Delta {\mu }_{\pm }}{\Delta \psi }\psi +e{z}_{\pm }\varphi ,$$

$${\vec{j}}_{\pm }^d$$

$${\vec{j}}_{\pm }^s$$

$${\vec{j}}_{\pm }^e$$

$${\vec{F}}_{I\pm }=\frac{k_{B}T}{D{c}_{\pm }}{\vec{J}}_{\pm }=-k_{B}T\nabla c_{\pm }-\frac{\Delta {\mu }_{\pm }}{\Delta \psi }{c}_{\pm }\nabla \psi -e{z}_{\pm }{c}_{\pm }\nabla \varphi .$$

$$\begin{array}{c}{\vec{F}}_{fs}=-\psi \nabla {\mu }_{\psi }=-\sum _{\pm }\psi \frac{\Delta {\mu }_{\pm }}{\Delta \psi }\nabla c_{\pm },\\ {\mu }_{\psi }=\frac{\delta {\mathcal{F}}_I}{\delta \psi }=\sum _{\pm }{c}_{\pm }\frac{\Delta {\mu }_{\pm }}{\Delta \psi }.\end{array}$$

$$\begin{split}\vec{F}_{f}=\vec{F}_{I+}+\vec{F}_{I-}+\vec{F}_{fs}=\\ -\sum_{\pm}\big{(}\underbrace{k_{B}T\nabla c_{\pm}}_{\vec{F}_{d}}+\underbrace{\frac{\Delta\mu_{\pm}}{\Delta\psi}\nabla(c_{\pm}\psi)}_{\vec{F}_{s}}+\underbrace{ez_{\pm}c_{\pm}\nabla\phi}_{\vec{F}_{e}}\big{)}.\end{split}$$

$$p=\frac{x_{0}^{2}}{3t_{0}^{2}}(n_{2}+n_{1})+\frac{x_{0}^{2}}{3}G\Psi_{1}\Psi_{2}+\sum_{\pm}\bigg{(}k_{B}Tc_{\pm}+\frac{\Delta\mu_{\pm}}{\Delta\psi}c_{\pm}\psi\bigg{)}.$$

$${\nabla }^2\rho /e$$

$$e{\nabla }^2{c}_t$$

$$e{\nabla }^2{c}_1$$

$$\to e{c}_1$$

$$\begin{array}{c}\frac{{\rho }_n}{e}=\sum _{i=1,3,5,..}^{2n+1}\frac{M^i{\Gamma }_{n-(i-1)/2}}{i!}{\Omega }^i+\sum _{j=2,4,6,..}^{2n}\frac{M^j{\mathrm{P}}_{n-j/2}}{j!}{\Omega }^j+{\mathrm{P}}_n,\\ e{c}_n=\sum _{i=1,3,5,..}^{2n-1}\frac{M^i{\mathrm{P}}_{n-(i-1)/2}}{i!}{\Omega }^i+\sum _{j=2,4,6,..}^{2n}\frac{M^j{\Gamma }_{n-j/2}}{j!}{\Omega }^j+{\Gamma }_n,\\ {\mathrm{P}}_n=-⟨{\rho }_n/e-{\mathrm{P}}_n⟩,{\Gamma }_n=-⟨c_n-{\Gamma }_n⟩+e{c}_s.\end{array}$$

$$\begin{gathered}\hat{c}_{1}=c_{s}\frac{\frac{M^{2}c_{s}}{2}\Omega^{2}+M\mathrm{P}_{0}\Omega}{\Upsilon_{1}},\\ \Upsilon_{1}=\bigg{\langle}\frac{M^{2}c_{s}}{2}\Omega^{2}+M\mathrm{P}_{0}\Omega\bigg{\rangle},\end{gathered}$$

$$\hat{\rho}_{1}/e=\frac{M^{3}c_{s}^{2}}{24\Upsilon_{1}}\Omega^{3}+\frac{M^{2}\mathrm{P}_{0}c_{s}}{2\Upsilon_{1}}\Omega^{2}\underbrace{-\bigg{\langle}\frac{M^{3}c_{s}^{2}}{24\Upsilon_{1}}\Omega^{3}+\frac{M^{2}\mathrm{P}_{0}c_{s}}{2\Upsilon_{1}}\Omega^{2}\bigg{\rangle}}_{=\mathrm{P}_{1}}.$$

$$\vec{F}_{s}=-\frac{\Delta\mu}{\Delta\psi}\nabla\Big{(}(c_{+}-c_{-})\psi\Big{)}=-2c_{s}M\frac{\Delta\mu}{\Delta\psi}\psi\nabla\psi.$$

$${\vec{F}}_{sc}^{\sigma }=-G{\Psi }_{\sigma }\nabla {\Psi }_{\bar{\sigma }}.$$

$$\begin{split}\psi\nabla\psi=2\psi\frac{n_{2}\nabla n_{1}-n_{1}\nabla n_{2}}{(n_{1}+n_{2})^{2}}=\\ =-2(n_{1}+n_{2})^{-3}\Big{(}n_{2}\nabla n_{1}+n_{1}\nabla n_{2}\Big{)}.\end{split}$$

$${\vec{F}}_d=-k_{B}T\nabla c_t=-k_{B}T{M}^2{\Gamma }_0\psi \nabla \psi /e.$$

$$\begin{array}{r}{\mathcal{F}}_T={\mathcal{F}}_F(\psi )+\int \xi |\nabla \Phi {|}^2\mathrm{d}\mathbf{r},\\ -{\nabla }^2\Phi =\psi .\end{array}$$

$$\begin{gathered}\Lambda=\frac{\lambda_{I}^{3}}{\epsilon\gamma}\bigg{(}\frac{ec_{s}\Delta\mu}{k_{B}T}\bigg{)}^{2},\\ \lambda_{ws}=2\pi\lambda_{I}\bigg{(}\frac{6}{\Lambda}\bigg{)}^{\frac{1}{4}}=2\pi\sqrt{\frac{k_{B}T}{ec_{s}\Delta\mu}}\big{(}6\lambda_{I}\epsilon\gamma\big{)}^{\frac{1}{4}},\\ \lambda_{ss}=4\lambda_{I}\bigg{(}\frac{12}{\Lambda}\bigg{)}^{\frac{1}{3}}=4\bigg{(}\frac{k_{B}T}{ec_{s}\Delta\mu}\bigg{)}^{\frac{2}{3}}\big{(}12\epsilon\gamma\big{)}^{\frac{1}{3}}.\end{gathered}$$

$$\gamma_{\rm{eff}}=\gamma-\int\xi|\nabla\Phi|^{2}\mathrm{d}z=\gamma\bigg{(}1-\frac{\Lambda}{2\lambda_{I}^{3}}\int\frac{|\nabla\Phi|^{2}}{\Delta\psi^{2}}\mathrm{d}z\bigg{)}.$$

$$\frac{|{\vec{j}}_{\pm }^s|}{|{\vec{j}}_{\pm }^e|}=4{\pi }^2\frac{{\lambda }_d^2}{{\lambda }_{ws}^2}\sqrt{\frac{W_i}{W_e}}=\frac{\Delta \mu }{\sqrt{6{\lambda }_I\gamma e^2/ϵ}}\sqrt{\frac{W_i}{W_e}}.$$

$$\begin{gathered}\Theta=\mathrm{arctan}\,\bigg{(}\frac{\partial_{x}\psi}{\partial_{y}\psi}\bigg{)}+\frac{\pi}{2}\\ \nu=\frac{\pi}{2}\frac{1}{\langle|\nabla\Theta|\rangle}.\end{gathered}$$

$$\begin{array}{c}{\lambda }_o=\frac{2\pi {\lambda }_d}{\sqrt{{\gamma }_p-1}}=\sqrt{\frac{4{\pi }^2ϵk_{B}T}{e^2{c}_s({\gamma }_p-1)}},\\ {\gamma }_p=\frac{|\frac{\Delta {\mu }_{+}}{k_{B}T\Delta \psi }-\frac{\Delta {\mu }_{-}}{k_{B}T\Delta \psi }|}{2\sqrt{{\kappa }_{sc}{e}^2/ϵ(k_{B}T{)}^2}}=\frac{\Delta \mu }{\sqrt{6{\gamma }_{sc}{\lambda }_I{e}^2/ϵ}},\end{array}$$

$$q_m=\frac{{\sum }_q{S}_{\psi }(q)q}{{\sum }_q{S}_{\psi }(q)}.$$

$$S_{\psi }^{-1}(q)\propto q^2+C{q}^{-2},$$

$$S_{\psi }^{-1}(q)\propto q^2+{\gamma }_p^2\frac{{\lambda }_d^{-4}}{q^2+{\lambda }_d^{-2}}.$$

Full text

††preprint: AIP/123-QED

I Introduction

Antagonistic salts are defined by the special property that, in a mixture of high and low permittivity fluids, one constituent ion kind migrates towards regions of higher permittivity while the other ion type does the opposite. Recent experimental studies of the antagonistic salt NaBPh4, solvated in a mixture of heavy water and 3-methylpyridine oil (3MP), suggest that adding antagonistic salt can cause structures to form on the nanometer scale Sadakane et al. (2013). These structures seem to be periodically repeating regions of high water/Na+ concentration and high oil/BPh${}_{4}^{-}$ concentration. Their length scale depends on many factors, such as the salt concentration Sadakane et al. (2013), volume fraction, and temperature Sadakane et al. (2009). Such mixtures show visible coloration with strong temperature sensitivity when the structure periodicity is of the order of the wavelength of visible light Sadakane et al. (2007). In other words, any fluid mixture of unequal permittivities may in principle be turned into a liquid crystal by simply adding an antagonistic salt. Being able to control mesoscopic structure formation in liquid crystals has potential for optical and nanoscale manufacturing technologies, such as directed self-assembly Müller and Rey (2018).

The theory of antagonistic salts has been studied to a notable extent by the group of Onuki Onuki et al. (2016). Their model describes the onset of structure formation in a system near the critical temperature of fluid demixing and the characteristic length of these structures Onuki and Kitamura (2004). Antagonistic salt mixtures have also been studied numerically by Araki, Okamoto and Onuki Araki and Onuki (2009); Onuki et al. (2011), showing lamellar structure formation and various degrees of nematic ordering. Recently, Tasios et al. Tasios et al. (2017) employed a simple lattice-based Monte Carlo algorithm to obtain a rich phase diagram of 3D structures from antagonistic salt mixtures. In the present work we add to this body of research an explicit analytical description of the coupling between fluid and ion concentrations and use it to show that antagonistic salt mixtures are very well approximated by the same free energy model commonly used to treat diblock copolymers – the Ohta-Kawasaki model Ohta and Kawasaki (1986). The Ohta-Kawasaki model is known to lead to ordered liquid crystal structures Zhang (2006); Weith et al. (2013); Cheng et al. (2017) seemingly identical to the structures we obtain numerically from antagonistic salt systems. By approximating our system with the Ohta-Kawasaki model we can predict for which parameters there is strong coupling of fluid and ion concentrations, which leads to structure formation, and when these structures become unstable due to zero or negative surface tension resulting from electrostatic contributions. We also add a more quantitative treatment of the dynamics of nematic ordering, showing that nematically ordered domains grow continuously at a rate proportional to the salt concentration.

This article is structured as follows. Section II explains our numerical method and justifies the simulation parameters used. We introduce a free energy model to derive the forces of fluid-ion coupling and demonstrate the relationship of our system to the Ohta-Kawasaki model. Our main results are split into five subsections. In Subsection III.1 we derive a polynomial function of the order parameter to approximate the charge distribution. Using this result we demonstrate in Subsection III.2 how our system is related to the Ohta-Kawasaki model. In Subsection III.3 we present our simulation results on nematic ordering. In Subsection III.4 we discuss and interpret experimental data from small-angle neutron scattering in light of our results and validate our model by comparing theoretical predictions of the characteristic length with the simulation results. Subsection III.5 concludes with a few examples of the different morphologies possible in 3D depending on volume fractions. Finally, in Sec. IV we summarize our results and discuss possibilities for future research.

II Method

We use the lattice-Boltzmann method (LBM) in an implementation previously presented by Rivas et al. Rivas et al. (2018); Rivas and Harting (2018). The LBM is an efficiently parallelizable numerical method to model fluid dynamics which is equivalent to the Navier-Stokes equation whenever flow velocities are much smaller than the speed of sound. In the LBM, space is discretized in a grid of lattice sites $x_{0}$ apart, each containing a set of populations $f_{i}.$ The origins of the LBM in statistical physics lead to the interpretation of each population as proportional to the probability of finding a fluid particle travelling in the direction given by the index $i$, but the fluid is modeled purely as a continuum nonetheless Benzi et al. (1992). The number of populations per lattice site is determined by the coarseness of the velocity discretization. We use 19 populations per lattice site (commonly called the D3Q19 scheme), corresponding to a zero vector representing resting particles plus 6 velocity vectors connecting the lattice site to its closest neighbours plus 12 to the next further neighbours in the grid Qian et al. (1992). The populations are evolved in each time step by an advection step, in which populations are simply moved to the lattice site they are connected with by their corresponding velocity vector, and by a collision step, in which viscous energy dissipation and fluid forcing are accounted for by relaxing the populations to an equilibrium distribution. In the single-relaxation-time scheme by Bhatnagar, Gross and Krook Bhatnagar et al. (1954) (BGK) that we use, the collision step takes the following shape:

[TABLE]

$f_{i}^{\mathrm{eq}}$ is a discretized Maxwell-Boltzmann distribution expanded in terms of Hermite polynomials up to second order in velocity Krüger et al. (2017); Benzi et al. (1992). $\tau$ gives the fraction of a time step, over which full relaxation to equilibrium takes place and determines the viscosity of the fluid. $S_{i}$, called the source term, contains the modification of the local equilibrium state due to any fluid forces added to the system. Our chosen method to calculate the source term goes back to Guo et al. Guo et al. (2002). A mixture of two fluids is simulated by simply representing each fluid component $\sigma$ by a separate set of populations $f^{\sigma}_{i}.$ The local mass concentrations of our fluid components are then $n_{\sigma}=\sum_{i}f^{\sigma}_{i}.$ We treat both fluids equally, i.e. they have equal density and viscosity.

The order parameter $\psi=(n_{1}-n_{2})/(n_{1}+n_{2})\in[-1,1]$ encodes the local fluid composition. We assume fluid incompressibility, i.e. $n_{1}+n_{2}=\mathrm{const.}=n_{\Sigma}.$ Fluid demixing is modeled by the pseudopotential method introduced by Shan and Chen Shan and Chen (1993). In this method, an interaction parameter $G$ encodes the strength of repulsive forces between unlike fluid components. These repulsive forces are included in the respective fluid’s source term $S^{\sigma}_{i}$ like any other force. For sufficiently high $G$, spinodal decomposition into two bulk phases occurs, whereas for sub-critical $G$ the phases mix homogenously. $G$ encodes both material properties and the effect of temperature on intermixability.

Experimental evidence of D2O/3MP mixtures shows structure formation induced by antagonistic salts at temperatures below the critical point, where macroscopic demixing would not occur without salt Sadakane et al. (2007). According to Onuki and Kitamura, this may be the effect of large thermal concentration fluctuations near the critical point of demixing stabilizing under the influence of the antagonistic salt Onuki and Kitamura (2004). As our model does not include thermal fluctuations, we only study the effect of antagonistic salt at values of $G$ slightly above the critical point of demixing.

The electrokinetic forces acting on fluids and ions are derived from a free energy functional $\mathcal{F}_{T}$ including fluid demixing via a phenomenological Ginzburg-Landau term $\mathcal{F}_{F}$ and ionic contributions $\mathcal{F}_{I}$ modelling the salt as an ideal gas coupled to the fluid by a solvation potential with an electrostatic contribution from the ionic charges.

[TABLE]

Here, $c_{\pm}$ is the concentration of positive/negative ions, $v_{0}$ is the ionic volume (assumed equal for both ion types for simplicity), $\rho=e(z_{+}c_{+}+z_{-}c_{-})$ is the charge concentration with the elementary charge $e$ and valence $z_{\pm}$, and $\phi$ is the electrostatic potential satisfying $\nabla^{2}\phi=-\rho/\epsilon.$ For simplicity, we assume the permittivity $\epsilon$ to be the same for both fluids and hence homogenous in space. We use a linearized solvation potential $\mu_{\pm}^{\mathrm{sol}}=\Delta\mu_{\pm}\psi/\Delta\psi$ with an ion-specific constant antagonicity $\Delta\mu_{\pm}$. $\Delta\mu_{\pm}$ is the solvation energy associated with transfering an ion from the bulk of one fluid phase to the other in a macroscopically demixed state where in the $n_{1}$-dominant phase $\psi=\psi_{1}$ and in the $n_{2}$-dominant phase $\psi=\psi_{2}$ fulfilling $\psi_{1}-\psi_{2}=\Delta\psi$. Note that throughout this work we write the Laplace operator as $\nabla^{2}$, while $\Delta$ signifies a variable related to a finite difference.

Similar models have been used in the past to theoretically describe binary fluid mixtures with antagonistic salts Onuki (2006); Rotenberg et al. (2009); Rivas et al. (2018). The fluid demixing term $\mathcal{F}_{F}=\int g(\psi,\nabla\psi)\hskip 2.15277pt\mathrm{d}V$ incorporates the surface tension $\gamma_{sc}$ caused by the repulsive pseudopotential forces in the LBM implementation. In order to model bulk demixing it should have the shape of a polynomial with two minima at the bulk values $\pm\Delta\psi/2$ of $\psi$ plus a gradient term. We choose $g(\psi,\nabla\psi)=g_{0}(\psi^{2}-\Delta\psi^{2}/4)^{2}+0.5\kappa_{sc}|\nabla\psi|^{2}.$ $g_{0}$ determines the energy penalty for mixed states. Its magnitude is irrelevant in our considerations. Minimization of $\mathcal{F}_{F}$ alone results in two bulk phases with an interface in the shape of a hyperbolic tangent $\psi(x)=0.5\Delta\psi\,\mathrm{tanh}\,(x/\lambda_{I}).$ We derive the interfacial width $\lambda_{I}=\Delta\psi^{-1}\sqrt{2\kappa_{sc}/g_{0}}$ by demanding that the chemical potential $\delta\mathcal{F}_{F}/\delta\psi=0$ at equilibrium. The surface tension is estimated by integrating the interfacial energy density from bulk to bulk, i.e. $\gamma_{sc}=0.5\int_{-\infty}^{+\infty}\kappa_{sc}|\nabla\psi|^{2}\,\mathrm{d}x=\Delta\psi^{3}\sqrt{\kappa_{sc}g_{0}/72}$, for a hyperbolic tangent shaped interface.

As long as spurious currents and inertial hydrodynamics are negligible, the pseudopotential model of demixing has been shown by Sbragaglia, Scarbolo et al. to correspond to a similar monotonically decreasing free energy functional Sbragaglia et al. (2009); Scarbolo et al. (2013). While full equivalence of the pseudopotential model to a free energy functional requires an additional gradient-shaped forcing term as a function of the order parameter, this term has been shown by numerical investigation to be generally negligible Sbragaglia et al. (2009).

The ion concentrations $c_{\pm}$ live on the same discrete lattice as the fluid and are evolved via a finite-difference scheme Rivas et al. (2018). By using a mean-field approach to model ion dynamics we neglect ion pair interactions, which may become important at salt concentrations high enough to cause significant steric interactions. The simplest possible time evolution of $c_{\pm}$ to conserve ion numbers and minimize the free energy in Eq. (2) is given by the Cahn-Hilliard equation Penrose and Fife (1990). In addition to the fluxes given by the Cahn-Hilliard equation we evolve the ions by an advection term in order to couple them to the velocity field $\vec{u}$ of the fluid mixture.

[TABLE]

$\mathcal{M}$ is the ion mobility, which can in general depend on the ion concentration and ion species. In order to recover Fick’s laws of diffusion we set $\mathcal{M}=D/k_{B}T$ according to the Einstein-Smoluchowski relation, with an ion diffusivity $D$ assumed to be constant. By performing a functional derivative of the ionic contributions $\mathcal{F}_{I}$ to the free energy, we obtain an expression for the ionic chemical potential:

[TABLE]

where we use $\phi=e/4\pi\epsilon\int(z_{+}c_{+}-z_{-}c_{-})/(|\vec{r}-\vec{r}_{0}|)\hskip 2.15277pt\mathrm{d}V$ to calculate the functional derivative of the electrostatic term. In solving Poisson’s equation, as in general, we treat our system as 3D, but with a thickness of one lattice site and periodic boundary conditions in the third dimension when emulating a 2D system. Using Eq. (4) we can write out the total ion flux $\vec{J}_{\pm}$ from Eq. (3). The resulting time evolution is identical to the well-known Nernst-Planck equation supplemented by a solvation term. The fluxes can be subdivided into diffusive, solvation, and electrostatic contributions $\vec{J}_{\pm}=\vec{j}^{d}_{\pm}+\vec{j}^{s}_{\pm}+\vec{j}^{e}_{\pm}$:

[TABLE]

Because the fluid is friction-coupled to the ions, it experiences the same force density as the ions. Assuming the inertial time scale of the ions to be much smaller than that of the fluid we expect the fluxes to correspond to the instantaneously reached drift velocity of the ions times the ion concentration and derive the force density experienced by the ions from Stokes’ law as

[TABLE]

In local equilibrium thermodynamics, spatial variations of intensive properties such as the chemical potential drive thermodynamic forces. The change in Gibbs free energy associated with a spatial variation of the chemical potential of fluid composition $\mu_{\psi}$ is $\mathrm{d}G=\psi\nabla\mu_{\psi}\mathrm{d}\vec{r}$, so that we can identify the corresponding force performing work on the system over an infinitesimal displacement $\mathrm{d}\vec{r}$ as

[TABLE]

The second term from the right in Eq. (6) proportional to $\nabla\psi$ represents the migration of ions towards regions of high concentration of the fluid species they are preferably solvated by. The fluid experiences this ion solvation force due to the aforementioned friction-coupling to the ions. Additionally the fluid experiences an analogous fluid solvation force representing migration of the fluid towards regions of high concentration of the ion kind preferably solvated by the local fluid composition. This force is not friction-coupled back to the ions because of the unequal inertial time scales of fluid and ions. The ions indirectly experience this force via the advection term in Eq. (3). In summary, the fluid forcing takes the form

[TABLE]

The equation of state, taking into account the gradient-shaped ideal ion pressure and total solvation forcing terms $\vec{F}_{d}$ and $\vec{F}_{s}$, but not the nonconservative electrostatic force $\vec{F}_{e}$, is then

[TABLE]

We write lattice units with subscript [math], i.e. the lattice length as $x_{0}$, one time step as $t_{0}$ and the units of mass as $m_{0}.$ In Eq. (9) we introduced the pseudopotential $\Psi_{\sigma}=m_{0}x_{0}^{-3}(1-\mathrm{exp}(-x_{0}^{3}n_{\sigma}/m_{0})).$ It approximates the fluid concentration of component $\sigma$, but is limited to values $\leq m_{0}x_{0}^{-3}$ (1 in simulation units). It is used in place of the fluid concentrations $n_{\sigma}$ only in order to avoid numerical instabilities, where fluid compression might increase the pseudopotential forces, in turn causing further fluid compression. The factors of 3 in Eq. (9) stem from the speed of sound of $1/\sqrt{3}x_{0}t_{0}^{-1}$ dictated by the chosen D3Q19 scheme of spatial discretization.

We restrict ourselves to monovalent ions $z_{+}=-z_{-}=1$ and symmetric antagonism $\Delta\mu_{+}=-\Delta\mu_{-}=\Delta\mu$ from here on. Simulation parameters roughly corresponding to experimental values can be chosen in the following way. First, we choose a length scale suitable to resolve the expected structure lengths, e.g. $x_{0}\approx 1$nm. The unit of mass is chosen so that $m_{0}x_{0}^{-3}$ matches the density of water. Due to a corresponding rescaling of the permittivity $\epsilon$ we can, without loss of generality, set the lattice unit of charge to the elementary charge $e_{0}=e.$

One limitation of the simple pseudopotential method we use is that the parameter $G$ controls both the interface width $\lambda_{I}$ and the surface tension $\gamma_{sc}.$ If we choose the time scale $t_{0}$ in such a way that the simulated fluid viscosity matches the viscosity of water, then for values of $\tau=1$, which are preferable in single-relaxation-time BGK for reasons of accuracy, we get a time step of $t_{0}\approx 0.2$ps. At such a small time step the surface tension obtained e.g. for $G=5.5x_{0}^{3}m_{0}^{-1}t_{0}^{-2}$ corresponds to 36 times the surface tension of water at room temperature. While one could lower $G$ to match the surface tension somewhat better, this would result in slower demixing dynamics, increased computational cost and extremely large interface widths of tens of nanometers. Instead, we choose $t_{0}$ in such a way, that we match the desired surface tension. As a side effect the viscosity in our simulations is significantly smaller than that of water, however we do not expect this to make any qualitative difference. The main source of fluid flow in our system is demixing, which is a diffusive process independent of inertial effects as long as $G$ is not chosen too large. While the resulting interface width is still of the order of a few nanometers, this is not neccessarily unrealistic for systems close to the critical point of demixing Buhn et al. (2004); Pousaneh et al. (2016), where structure formation by antagonistic salts is observed in experiments. When approaching the critical point, the interface length $\lambda_{I}$ diverges and the surface tension $\gamma_{sc}$ vanishes with a strong temperature sensitivity.

Choosing for example $x_{0}=$1nm, $t_{0}=10$ps and $G=5.5x_{0}^{3}m_{0}^{-1}t_{0}^{-2}$ we have $\gamma_{sc}=0.036m_{0}t_{0}^{-2}=2.75\cdot 10^{-4}$N m*-1* and $\lambda_{I}=1.3x_{0}=1.3$nm. Accordingly, $\epsilon_{r}=40$ corresponds to $\epsilon=0.1$ in simulation units, $k_{B}T=0.6m_{0}x_{0}^{2}t_{0}^{-2}$ for a temperature of about 330K and a salt concentration of 0.1 is equivalent to about 166mmol l${}^{-1}.$ In practice, as our aim is to understand the system fundamentally, we experimented with a wide range of simulation parameters without focusing very much on their relation to experimental values. For example in most simulations we keep $\Delta\mu$ in a range of about $3-8k_{B}T$, which is noticeably lower than the $15k_{B}T$ of NaBPh4 commonly used in experiments, because this allows us to more easily avoid numerical errors from steep ion concentration gradients and high fluid forcing terms. Doing so we have come to the conclusion that, above all, the dimensionless numbers $\Lambda$ and $\lambda_{d}/\lambda_{ws}$ introduced in Eqs. (21) and (23) determine the behaviour of the system. In the interest of reproducibility, we nonetheless give the simulation parameters used in simulation units in all figure captions.

III Results

III.1 Fluid-charge coupling

All our simulations are initialized as a homogenous mixture with symmetric volume fraction $m=\langle n_{1}\rangle/(\langle n_{1}\rangle+\langle n_{2}\rangle)=0.5$ and constant ion concentrations $c_{\pm}=0.5c_{s}$ everywhere given by the initial salt concentration $c_{s}.$ The fluid concentrations $n_{1}$ and $n_{2}$ are initially perturbed by a small stochastic prefactor at each lattice site in a way that keeps the overall volume fraction constant. For $G=5.5x_{0}^{3}m_{0}^{-1}t_{0}^{-2}$ spinodal decomposition during the first about one thousand time steps is driven almost completely by the pseudopotential forces of demixing. As the local contrast in fluid concentrations $n_{1}$ and $n_{2}$ grows, charges build up due to ion separation by solvation and the resulting electrostatic forces on the fluid arrest the process of coarsening. After a further about one thousand time steps we observe lamellar, bicontinuous structures.

Because of a strong coupling via the solvation potential, both the fluid order parameter $\psi$ and the charge $\rho$ follow the same lamellar pattern, as shown in Fig. 1. Our key result in this paper is an analytic description of this coupling. Using the linear coupling limit valid for low $\Delta\mu/k_{B}T$ allows us to use prior results on the Ohta-Kawasaki free energy model to derive the two dimensionless numbers $\Lambda$ and $\lambda_{d}/\lambda_{ws}$ describing our system.

As we are concerned with systems close to the critical point of demixing, the pseudopotential forces of demixing are small and the time scale of fluid relaxation is much larger than that of ion relaxation. This leaves the advection term in Eq. (3) largely irrelevant. In our simulations the ions are practically in equilibrium at all times with respect to the current fluid composition $\psi.$ From the equilibrium solution of the evolution Eq. (3) for the charge $\rho$ and the total local ion concentration $c_{t}=c_{+}+c_{-}$, we obtain two coupled differential equations,

[TABLE]

For ease of notation we write $M=-\Delta\mu/\Delta\psi k_{B}T$ and $\Omega=\psi+e\phi\Delta\psi/\Delta\mu.$ When $M\lesssim 1$, solvation is weak compared to the ideal pressure of ions and diffusion keeps the total ion concentration homogenous, so that $c_{t}\approx c_{0}=\mathrm{const.}=c_{s}=\Gamma_{0}.$ Using $c_{t}=c_{0}$ our first approximation of the charge follows from Eq. (10) as $\rho_{0}/e=c_{s}M\Omega-\langle c_{s}M\Omega\rangle=c_{s}M\Omega+\mathrm{P}_{0}$, where $\mathrm{P}_{0}$ is chosen to ensure charge neutrality. Here, as throughout the rest of the paper, $\langle\cdot\rangle$ signifies a spatial average. $\rho_{0}$ is a valid approximation as long as $\rho$ does not depend linearly on spatial coordinates, which might be the case, e.g., if $\phi$ includes strong external electric fields. By inserting this into Eq. (11), we approximate the first order fluctuations of $c_{t}$ around its expectation value $c_{s}$, which will be significant for larger $M$:

[TABLE]

Here, $\Gamma_{1}$ is chosen to ensure total ion conservation. The approximation $c_{t}=c_{1}$ can in turn be plugged back into Eq. (10), yielding a higher-order approximation $\rho_{1}$ for the charge $\rho$, and so on. The general expansions of $\rho$ and $c_{t}$ up to order $n$ are

[TABLE]

It is worth noting that $\langle\Omega^{k}\rangle=0$ for uneven $k$ whenever the volume fraction is $1/2$ and the two fluids are treated identically apart from $\Delta\mu_{+}=-\Delta\mu_{-}$ and $z_{+}=-z_{-}.$ This in turn means that all $P_{n}=0$ and thus Eq. (13) becomes much simpler, with $\rho_{n}$ turning into a polynomial of only uneven powers of $\Omega$ and $c_{n}$ of only even powers of $\Omega.$

Eq. (13) fails by overestimating the fluctuations of $c_{t}$ when $\Delta\mu\gg k_{B}T.$ Strong solvation leads to completely ionless regions of $c_{t}\approx 0$ in the vicinity of the interfaces, because here both ion kinds experience solvation forces in opposite directions towards higher concentrations of their preferred fluid components. Numerically, whenever outgoing fluxes at any lattice site would reduce the local ion concentration below zero, all outgoing fluxes at that site are reduced by a common factor, so that the local ion concentration goes to zero instead. Ingoing fluxes are then calculated and potentially similarly reduced individually by the previously computed factor reducing the neighbouring lattice site’s outgoing fluxes. In our simulations, ion concentrations may become negative without this type of discretization correction when $\Delta\mu\gtrsim 7k_{B}T.$

The analytic approximations leading to Eq. (13) do not take a limited ion availability into account and similarly predict negative $c_{t}$ for $\Delta\mu\gg k_{B}T.$ To a surprisingly good degree of accuracy we can still approximate $c_{t}$ in this limit by rescaling each consecutive $c_{n}$ by a factor $\Upsilon_{n}$ in order to force an amplitude of ion concentration fluctuations equal to $c_{s}.$ The additive constants $\Gamma_{n}$ then all become [math] for $n\geq 1.$ We write the charge and ion concentration approximations including the correction for finite salt concentration as $\hat{\rho}_{n}$ and $\hat{c}_{n}.$ Thus Eq. (12) becomes

[TABLE]

and by reinserting into Eq. (10) we obtain the charge as before.

[TABLE]

Naturally, the procedure can be repeated to obtain higher order corrections and a general expansion as in Eq. (13). While $\rho_{n}$ and $c_{n}$ tend to converge towards the exact results $\rho$ and $c_{t}$ for $n\to\infty$ as long as $\Delta\mu$ is small enough that ion-depletion at the interfaces does not play a role (roughly $\Delta\mu<5k_{B}T$), the ideal value of $n$ for $\hat{\rho}_{n}$ and $\hat{c}_{n}$ must be empirically chosen depending on $\Delta\mu$. We find that $n=2$ produces a more accurate approximation than $n=1$ only for high values of $\Delta\mu\approx 15k_{B}T$, the largest antagonicity we used in our simulations. Higher values of $n$ would likely only be appropriate for extreme antagonicities of $\Delta\mu\gg 15k_{B}T.$

The proposed methods of approximating the ion distributions gain their value from the possibility of neglecting electrostatic fluxes, so that $\Omega\approx\psi.$ In Fig. 2 we show the quality of the theoretical model using $\Omega=\psi$ for a range of values of $c_{s}$ and $\Delta\mu.$ In Fig. 2a the root mean square errors of approximating $\rho$ as $\rho_{n}$ are seen to decrease with increasing $n$ and $c_{s}.$ In Fig. 2c we can see that, for low $c_{s}$ and thus large structure sizes, the charge begins to concentrate at the interfaces due to electrostatic interaction with the counter charges in the opposite fluid component. This feature is not captured by $\rho_{n}$ due to neglecting electrostatics.

Errors depend in a more complicated manner on $\Delta\mu$ than on $c_{s}.$ As shown in Fig. 2b, $\rho_{2}$, by ignoring limited availability of ions, quickly fails for high $\Delta\mu$ by predicting negative ion concentrations. For $\Delta\mu\to 0$ it fails again because structure sizes diverge when there is no solvation, causing $\rho$ and $\psi$ to become decoupled. The charge approximations $\hat{\rho}_{n}$ corrected for limited ion availability fare much better for high $\Delta\mu$, but the errors here do not converge to 0 for $n\to\infty.$ Instead, $\hat{\rho}_{1}$ is a better approximation for intermediate $\Delta\mu$ and $\hat{\rho}_{2}$ for high $\Delta\mu.$ In Figs. 2e and f we can see how for larger $n$, $\hat{\rho}_{n}$ features more pronounced regions of zero charge at the interfaces, which correspond to ion-depleted regions of $c_{t}\approx 0.$

III.2 Ohta-Kawasaki equivalency

Using the analytical results for the charge and ion concentrations from Sec. III.1 we can significantly simplify our system and make use of known results from the Ohta-Kawasaki free energy model. Looking back at the solvation force in Eq. (8) we find, by taking the linear approximation for symmetric fluids $\rho=\rho_{0}=ec_{s}M\psi$

[TABLE]

This is very similar to the continuum limit of the pseudopotential force acting on fluid $\sigma$ due to interaction with fluid $\bar{\sigma}$ in the pseudopotential model when neglecting third and higher order gradient terms Sbragaglia et al. (2009):

[TABLE]

Setting $\Psi_{\sigma}=n_{\sigma}$ and postulating $\nabla n_{2}=-\nabla n_{1}$ due to incompressibility, we find that $\vec{F}_{s}$ in fact points in precisely the same direction as $\vec{F}^{1}_{sc}+\vec{F}^{2}_{sc}$, that is, away from the interface.

[TABLE]

This means that the solvation force acts purely to increase surface tension. When fluid incompressibility is fulfilled, $n_{1}+n_{2}$ is constant and Eq. (18) has the same shape as Eq. (17). Similarly, the ideal pressure force resulting from ion diffusion also acts only in the direction of the pseudopotential force assuming symmetric fluids and therewith $c_{t}=c_{1}=\frac{M^{2}\Gamma_{0}}{2e}\psi^{2}+\Gamma_{1}/e.$

[TABLE]

Indeed we find both $\vec{F}_{d}$ and $\vec{F}_{s}$ to be almost irrelevant in determining the general morphology of the system in simulations as long as their magnitude does not rise to a level where they cause fluid compression. This holds true even when $\vec{F}_{d}$ and $\vec{F}_{s}$ are orders of magnitude larger than the electrostatic force $\vec{F}_{e}.$ By subsuming the influence of $\vec{F}_{d}$ and $\vec{F}_{s}$ into a modified surface tension $\gamma=\gamma_{sc}+\gamma_{d}+\gamma_{s}$ and neglecting the electrostatic fluxes $\vec{j}_{\pm}^{e}$, we are left with a greatly simplified system. The surface tension $\gamma$ can be determined from Eq. (9), e.g. via a Laplace test, meaning that a single bubble of one fluid component of various radii is initialized in the bulk of the other fluid component. According to the Young-Laplace law, the surface tension $\gamma$ is the slope of the graph of pressure differences inside to outside the droplet versus the inverse droplet radius.

Whenever $\Delta\mu$ is small enough for the first-order charge approximation $\rho_{0}/e=c_{s}M\psi$ to be accurate, our system is fully equivalent to the Ohta-Kawasaki model commonly used to model diblock copolymers Spadaro (2009):

[TABLE]

The first term is the fluid demixing energy $\mathcal{F}_{F}$ introduced in section II but with $\kappa_{sc}$ substituted with $\kappa$, which fulfills $\gamma=0.5\int_{-\infty}^{+\infty}\kappa|\nabla\psi|^{2}\,\mathrm{d}x$ for a hyperbolic tangent shaped $\psi$ and thus includes the surface tension contributions of $\vec{F}_{d}$ and $\vec{F}_{s}$. Crucially the model adds an electrostatic term weighted by $\xi$ in which the order parameter directly corresponds to a charge. The long-range interaction potential $\Phi$ is, in the first-order charge approximation $\rho=\rho_{0}$, proportional to the electrostatic potential with $\Phi=\phi\epsilon/ec_{s}M.$

The Ohta-Kawasaki model is usually treated in either of two limits depending on the characteristic length $\lambda_{L}$ of structures in equilibrium as compared to the interface width $\lambda_{I}.$ In the weak-segregation limit, i.e. when $\lambda_{L}\lesssim\lambda_{I}$, the wavelength of structures is Weith et al. (2013) $\lambda_{ws}=2\pi(\kappa/2\xi)^{1/4}$, or, as $\kappa=6\gamma\lambda_{I}/\Delta\psi^{2}$, $\lambda_{ws}\propto(\gamma/\xi)^{\frac{1}{4}}.$ In the strong-segregation limit Weith et al. (2013), i.e. when $\lambda_{L}\gg\lambda_{I}$, $\lambda_{ss}\propto(\gamma/\xi)^{\frac{1}{3}}.$ In analogy to our free energy Eq. (2), the term $\xi|\nabla\Phi|^{2}$ corresponds to the electrostatic term $0.5\rho\phi.$ By setting $\rho=\rho_{0}=ec_{s}M\psi$ we can see by comparison that $\xi=(\Delta\mu/k_{B}T\Delta\psi)^{2}e^{2}c_{s}^{2}/2\epsilon=\Lambda\gamma/2\Delta\psi^{2}\lambda_{I}^{3}.$ In summary:

[TABLE]

The dimensionless number $\Lambda$ turns out to be useful to predict the morphology produced by a given set of system parameters. In our simulations we find by varying all involved simulation parameters, including the pseudopotential interaction parameter $G$, that structures begin to gradually dissolve in the range of $\Lambda\approx 1-10$, or $\lambda_{ws}/\lambda_{I}\approx 6-9$ up to a return to a mixed state at higher $\Lambda.$ The onset of structure dissolution is characterized by lamellar regions interspersed with blots of mixed regions.

Following Onuki and Okamoto Onuki and Okamoto (2009) the electrostatic surface tension contribution in a system with heterogenity in only the $z$ direction, such as parallel lamellar structures, can be calculated by subtracting the electrostatic energy density integrated from one bulk to the other over an interface:

[TABLE]

By making a simple 1D ansatz of $\psi=0.5\Delta\psi\,\mathrm{sin}(2\pi x/\lambda_{ws})$ valid in the weak-segregation limit and integrating $|\nabla\Phi|^{2}$ from the bulk of one phase to the other over $\lambda_{ws}/2$, we can predict a critical $\Lambda=2^{13}3^{-3}\pi^{-4}\approx 3.1$, at which $\gamma_{\rm{}_{eff}}$ is expected to become zero. This approximately matches the point at which we observe the onset of structure dissolution in simulations.

Recall that in deriving the equivalency to the Ohta-Kawasaki model we set $\Omega\approx\psi$, i.e. we neglected electrostatic fluxes. Defining the interfacial energy density $W_{i}=0.5\kappa|\nabla\psi|^{2}$ and the electrostatic energy density $W_{e}=0.5\epsilon|\nabla\phi|^{2}$ we can, using again $\rho=\rho_{0}$, quantify the ratio of solvation and electrostatic fluxes as

[TABLE]

Here we introduce the Debye length $\lambda_{d}=\sqrt{\epsilon k_{B}T/e^{2}c_{s}}.$ It gives the length scale over which the density of countercharges near a surface charge decays. The interfacial and electrostatic energy densities are the main two counteracting factors driving structure formation in the Ohta-Kawasaki model, with the former striving to minimize interface area and the latter favouring either homogenous charge, or rapid spatial oscillations of $\rho$ and therewith $\psi$, leading to a large number of interfaces. A balance of these two energy densities is found when lamellar structures form, so it seems reasonable to assume a comparable order of magnitude in such morphologies. Indeed, Araki and Onuki showed this to be true in steady-states and weak segregation Araki and Onuki (2009). Assuming $W_{i}\approx W_{e}$, Eq. (23) suggests electrostatic fluxes to be about an order of magnitude smaller than solvation fluxes and hence quantitatively negligible whenever the Debye length is $\gtrsim\lambda_{ws}/2.$ In other words, for $\lambda_{d}\gtrsim\lambda_{ws}/2$ ion dynamics are largely unaffected by electrostatic interactions with other ions. When $\lambda_{d}\to\lambda_{ws}/2$ structure sizes diverge and for $\lambda_{d}<\lambda_{ws}/2$ spinodal decomposition is largely unaffected by electrostatic effects.

III.3 Nematic ordering

As long as the ion dynamics are dominated by solvation, the morphology is almost completely determined by $\Lambda.$ As $\Lambda$ is increased, lamellae become thinner and the degree of nematic ordering at a given time increases. In the top row of Fig. 3 we show three representative morphologies for $\Lambda=0.08$, 0.9 and 1.3.

Nematic order parameters commonly used in the analysis of liquid crystal structures measure the overall homogenity of nematic orientation. They are suitable when one is interested in deviations from some prefered direction, for example given by an applied electric field, surface patterning or initially ordered state Buluy et al. (2018); Luders et al. (2015); Saintillan and Shelley (2007), or to determine the nematic homogenity of some spontaneously ordered steady-state Purdy et al. (2003); Ginelli et al. (2010). In our case, ordering progresses without any inherently preferred direction and we are interested in the time evolution of the size of nematically ordered domains. As long as these domains are much smaller than the simulation domain, standard nematic order parameters indicate no nematic order, because all nematic orientations are equally frequent in the simulation domain. When the size of nematically ordered domains approaches the simulation domain size, some random nematic orientation begins to dominate and standard nematic order parameters indicate order, but we are uninterested in this regime because under these conditions finite size effects affect the time evolution of nematic order.

Instead we quantify nematic ordering over time by a nematic range parameter $\nu$ giving the length scale over which the nematic orientation changes. For this purpose we first determine the angle of orientation $\Theta$ of the lamellae at a given position. The bottom row of Fig. 3 shows the local angle of orientation as a function of position for three morphologies. The gradient magnitude of the orientation field $\Theta$ gives the average angle by which the lamellar orientation changes over a distance of one lattice site in the direction of the greatest rate of change. The local nematic range $\nu$ estimates the distance in lattice sites in direction of fastest change of $\Theta$ over which $\Theta$ changes to a perpendicular orientation. Using the fact that the gradient of $\psi$ ought to be everywhere perpendicular to the local lamella, we choose the definition

[TABLE]

When the entire system domain is filled with exactly parallel lamellae, $\nabla\Theta$ is zero everywhere, and $\nu\to\infty.$ Otherwise, the smaller the domains of common orientation are, the larger $|\nabla\Theta|$ is on average, so that $\nu\to 0$ for disordered systems.

The driving force of nematic ordering is long-range electrostatic fluid forcing. As our model does not include thermal fluctuations, no opposing force exists to disrupt nematic ordering. Hence, nematic ordering appears to continue indefinitely, until either the entire system domain is filled by perfectly ordered lamellae, or the process of nematic ordering is hindered by finite size effects. In Figs. 4a and c, we show how $\nu$ increases mostly monotonously in time. Nematic ordering begins to fluctuate considerably when domains of common orientation become comparable in size to the simulation domain, which occurs in Fig. 4a when $\nu\gtrsim 38x_{0}.$ While this corresponds to only about 8% of the system size of $504x_{0}$, recall that $\nu$ measures the size of ordered domains in the direction of the fastest change of $\Theta.$ In any other direction, the domain may be considerably larger. One example of finite size effects contributing to fluctuations of $\nu$ is the orientation angle dependence of the number of separate lamellae that fit into the simulation domain due to its square shape and the use of periodic boundary conditions. For this reason large numbers of nematic defects must form, leading to decreasing $\nu$, before nematic ordering can increase further.

As nematic ordering is driven by electrostatic fluid forcing, and the electrostatic fluid forcing is proportional in magnitude to $\Lambda$, it is not surprising that ordering progresses faster for higher $\Lambda.$ Nematic range as a function of time seems to be well-modelled as a power law of the form $\nu/x_{0}=(t/t_{0})^{\eta}$, with an exponent $\eta\propto c_{s}$ proportional to the salt concentration. Because of the early onset of finite size effects we made another simulation on a much larger 2064$x_{0}$ system over $3\cdot 10^{6}$ time steps for $\Lambda=1.16$. The resulting evolution of $\nu$ shown in Fig. 4c confirms the previously mentioned power law behaviour but without any major deviations due to finite size effects. Such a power law behaviour has been previously described for the growth of nematically ordered domains in the Ohta-Kawasaki model Huang and Viñnals (2007). As shown in Eq. 21, $\Lambda\propto c_{s}^{2}.$ Judging from smaller-scale simulations, the system parameters $\epsilon$ and $\Delta\mu$ affecting the magnitude of electrostatic forces and therewith $\Lambda$ appear to similarly affect the rate of nematic ordering. We suspect that indeed $\eta\propto\sqrt{\Lambda}$, though we have so far only performed sufficiently long-term simulations to accurately determine $\eta$ as a function of the salt concentration. The results of fitting the exponents $\eta$ as either $\propto c_{s}$ or $\propto\sqrt{\Lambda}$ are shown in Fig. 4b. We are currently engaged in implementing a performance efficient version of our code for 2D simulations using the D2Q9 scheme in order to further study nematic ordering as a function of $\Lambda$, for larger system sizes and time scales, and to test the influence of viscosity. It is known that in the Ohta-Kawasaki model, $\eta$ is a function of $\xi\propto\Lambda$, though the precise form of this dependency has not been studied to our knowledge Huang and Viñnals (2007).

Because the lamellar patterns are oppositely charged, the nematic ordering can be significantly sped up and oriented in a controlled manner by applying an external electric field. The electric field has to be applied early during the demixing process to be able to align the lamellae. When an electric field is applied on a domain of already nematically ordered lamellae, the resulting electrostatic forces acting on both sides of an interface are directed in opposition to each other and cancel out due to opposite charges. Increasing the field strength eventually leads to a break-up of the lamellar structures which is known as the Helfrich-Hurault instability in the context of smectic and cholesteric liquid crystals Onuki and Fukuda (1995). In the case of our simulations, hydrodynamic chaos ensues, though in the Stokes regime, a square lattice undulating pattern can be predicted from free energy considerations Onuki and Fukuda (1995); Fukuda and Onuki (1995). If the field strength is decreased again, lamellar structures in alignment with the external field direction will reform. Structure formation by electric fields has been extensively studied for diblock-copolymer mixtures Böker et al. (2002); Kan and He (2016); Pester et al. (2017) with results that should be entirely transferable to the case of antagonistic salt mixtures as long as external electric fields are not strong enough to decouple the ions from the fluid order parameter.

III.4 Length scales

In small-angle neutron scattering (SANS) a neutron diffraction image of a sample is recorded and time-averaged on a screen. According to the Rayleigh-Gans equation, this diffraction image is well-approximated as the radially averaged power spectral density (PSD) of the scattering length density in the sample. The scattering length density is an empirically determined material constant. Fig. 5a shows the PSD $S_{n_{s}}(q)$, i.e. the absolute squared of the Fourier transform, of a linear combination $n_{s}=s_{1}n_{1}+s_{2}n_{2}+s_{+}c_{+}+s_{-}c_{-}$ of the fluid and ion concentrations. We choose the prefactors $s_{\sigma}$ and $s_{\pm}$ representing relative scattering length densities as $s_{1}=1.4$, $s_{2}=6.4$, $s_{-}=2.1$ and $s_{+}=0$ from the experimental values for a D2O/3MP mixture with the antagonistic salt NaBPH4 Sadakane et al. (2013). Note that we do not aim to quantitatively match the experimental data, as we are operating at different volume fractions and neglecting potentially important system parameters such as permittivity and viscosity differences between the fluids. Nonetheless, $S_{n_{s}}(q)$ has a notable resemblance to experimental results from SANS Sadakane et al. (2013). For low salt concentrations, or low $\Lambda$, the PSD approximately follows the Ornstein-Zernike function Ornstein and Zernike (1914); Sadakane et al. (2013), indicating typical hydrodynamic concentration fluctuations but no periodic structure. At intermediate salt concentrations, a major peak appears at a wavenumber $q_{m}=1/\lambda_{L}$ corresponding to the center-to-center spacing $\lambda_{L}$ between two neighbouring lamellae of either fluid.

The secondary peaks at higher wavenumbers for intermediate $\Lambda$ in Fig. 5a are the harmonics of the peak frequency, i.e. they are located at integer multiples of $q_{m}.$ The uneven harmonics can be explained by approximating the lamellar structures as a square wave of frequency $q_{m}$, which by the Fourier series can be decomposed into the sum of all uneven harmonics of $q_{m}.$ The even harmonics are caused by a decrease of the total fluid density at the interfaces due to the pseudopotential forces of demixing. Our general assumption of fluid incompressibility is of course only approximately true here. The density dip at the interface can also be approximated as a step function, but at twice the frequency $q_{m}$ of the lamellae, as each lamellar has two interfaces. At least one such secondary peak at quite precisely double the frequency of the first peak is visible in the experimental data of Sadakane et al. Sadakane et al. (2013), and as its intensity is low compared to the noise, we believe that higher order peaks may have simply not been resolved due to imprecisions of the measurement. The PSD of the order parameter $S_{\psi}(q)$, while otherwise almost identical to $S_{n_{s}}(q)$, does not show secondary peaks at the even harmonics of $q_{m}.$ The sharp interfacial density dips in $n_{1}$ and $n_{2}$ are smoothed by normalization with $n_{1}+n_{2}.$

Going to high salt concentrations, we again recover essentially the same picture in the PSD as for low salt concentrations. As we illustrate in Fig. 5b, periodic structures destabilize and remixing occurs, when $\Lambda$ exceeds some threshold. In section III.2 we estimated the critical value where the effective surface tension goes to zero as $\Lambda_{c}\approx 3.1.$ In simulations the point where we observe structure dissolution varies depending on the values of $G$ and $\Delta\mu/k_{B}T$ in the range of about $\Lambda_{c}\approx 1-10.$ A likely reason for this is that the sinusoidal single-wave approximation we made in estimating $\Lambda_{c}$ is overly simplistic. The actual shape of the lamellae can have aspects of a square wave as well as flat shoulders of almost constant zero-valued order parameter at the interfaces for high $\Delta\mu$, as shown in Figs. 2e and f. Also, an imperfect nematic ordering changes the electrostatic field in a very non-trivial way.

It is worth noting that the remixing state shown on the right in Fig. 5b does not ever reach fluid flow equilibrium, with domains of non-zero order parameter sprouting and dissolving continuously. We expect that this likely unphysical behaviour might disappear when using a less coarse-grained model, where each fluid molecule experiences only the electrostatic force acting on ions it is currently bound to by solvation. As it is, electrostatic forces are applied equally to both fluid components at each lattice site.

Extracting the structure size, i.e. the periodicity $\lambda_{L}$ of the lamellae as the inverse of the spatial frequency of the primary peak in the PSD we find good agreement with the scaling laws derived in Eq. (21) according to the Ohta-Kawasaki model. We compare the theoretical predictions with simulation results for varying salt concentrations, dielectric permittivity and antagonicity in Figs. 6a and b. The parameters $\lambda_{I}$ and $\gamma$ are determined via Laplace tests. $\lambda_{I}$ is determined by fitting a hyperbolic tangent function to a 1D cut through the droplet interface. Full lines show the predicted structure sizes using a value of $\gamma$ obtained from Laplace tests including ionic contributions. Here the surface tension is calculated using Eq. (9), so that ion pressure and solvation effects are present but electrostatic contributions are disabled by setting $\vec{j}^{e}_{\pm}=\vec{F}_{e}=0$ in the Laplace test. For dashed lines the Laplace test is performed without any ions present, so that $\gamma=\gamma_{sc}.$ Although $\gamma$ changes by about a factor of 2 due to ionic contributions from the lowest to the highest salt concentration, the impact on the structure size predictions is not very large. Non-electrostatic ionic contributions to surface tension may be neglected in calculating the structure size unless a particularly high degree of accuracy is desired.

Following Onuki and Kitamura Onuki and Kitamura (2004), the structure size is given by

[TABLE]

with a dimensionless parameter $\gamma_{p}$ quantifying the strength of antagonicity. Our weak-segregation scaling in Eq. (21) can be rewritten as $\lambda_{ws}=2\pi\lambda_{d}/\sqrt{\gamma_{p}}$, which is almost identical to Onuki’s prediction. The model by Onuki predicts a divergence of structure sizes when $\gamma_{p}\to 1$, which our model does not reproduce directly. The reason for this is that the Ohta-Kawasaki model, from which the scaling laws in Eq. (21) are derived, describes our system accurately only when electrostatic fluxes are small versus solvation fluxes and thus $\gamma_{p}\gg 1$ (cf. Eq. (23)). When $\gamma_{p}\lesssim 1$, electrostatic fluxes keep the ions bound to the interfaces, the charge is no longer strongly coupled to the order parameter and periodic structure formation ceases. We find Onuki’s prediction to be almost identical to the weak-segregation Ohta-Kawasaki limit in Figs. 6a and b for $\gamma_{p}\approx 11.$ Onuki’s model fares slightly worse than the strong-segregation Ohta-Kawasaki limit when structure sizes grow larger than about $15\lambda_{I},$ but Onuki’s model can be expected to fare better when $\gamma_{p}\to 1.$

The position of the primary spatial frequency peak, i.e. the structure size, can also be quite accurately determined by calculating the radially averaged PSD of the order parameter $S_{\psi}(q)$, i.e. the structure factor, and then taking its first moment:

[TABLE]

The dominant length scale taken as the inverse of $q_{m}$ converges very quickly compared to the slow convergence of nematic ordering, as seen in Fig. 6c. Note the slower convergence for lower salt concentrations, as it is the electrostatic forcing, scaling with the salt concentration, that stops the spinodal decomposition in the first place. As structure sizes are increased at a constant interface width $\lambda_{I}$ we gradually approach the limit of strong segregation, in which Araki and Onuki similarly observed a noticeably slower convergence of structure sizes Araki and Onuki (2009).

Comparing the structure factor of the order parameter in the model by Onuki and Kitamura Onuki and Kitamura (2004) with that in the model of Ohta and Kawasaki Ohta and Kawasaki (1986) reveals why a divergence of structure sizes occurs in the former as $\gamma_{p}\to 1$ but not in the latter. In the Ohta-Kawasaki model the inverse structure factor can be approximated as

[TABLE]

with some constant $C>0$ and a $q^{-2}$ term stemming from Coulombic attraction of opposite fluid phases inhibiting macroscopic demixing and forcing $S_{\psi}\to 0$ for $q\to 0$. In the description of pure diblock copolymer solutions this is of course a reasonable condition indicating that oppositely charged blocks of a single diblock copolymer cannot stretch and separate indefinitely Ohta and Kawasaki (1986); Schmid (2011). In the model of Onuki and Kitamura on the other hand Onuki and Kitamura (2004)

[TABLE]

Due to the presence of the $\lambda_{d}^{-2}$ term, we no longer neccessarily have $S_{\psi}\to 0$ for $q\to 0$. The fact that long-range interactions of charges are screened by the Debye layer allows for macroscopic demixing. By performing the first and second derivatives of $S_{\psi}^{-1}(q)|_{q=0}$ we find that the structure factor $S_{\psi}(q)$ has a maximum at $q=0$ for $\gamma_{p}\leq 1$ and a minimum at $q=0$ for $\gamma_{p}>1$. Thus macroscopic demixing occurs for $\gamma_{p}\leq 1$. A similar form of structure factor as in the model of Onuki and Kitamura can also be used to describe polyelectrolytes stabilized by electrostatic repulsion in a poor solvent. In such systems mesophase structures are formed for small salt concentrations and a transition to macroscopic demixing occurs when $\lambda_{d}$ falls below some threshold value Joanny and Leibler (1990); Raphael and Joanny (1990).

III.5 Extension to 3D

The method as discussed is extensible to 3D in a straightforward manner. Preliminary results so far are essentially identical to the 2D case, with gradual nematic ordering in proportion to the salt concentration, a charge distribution well-approximated as a polynomial function of the order parameter and average structure sizes as predicted by the Ohta-Kawasaki and Onuki models in Eq. (21). A somewhat wider variety of different periodic structures is observed depending on the volume fraction $m$, as shown in Fig. 7. A minority phase tends to form spherical bubbles, which may split or elongate depending on the strength of electrostatic forces. Slightly asymmetric volume fractions lead to tube structures, which may be separate or, for almost symmetric volume fractions, conjoined into a bicontinuous network. Lastly, symmetric volume fractions lead to lamellar structures. The same morphologies as a function of volume fraction have been previously produced by the Ohta-Kawasaki model Thomas (2007). In 2D we observe only lamellar and droplet phases.

An interesting avenue of further research lies in the transition region between the droplet and tubular phases. Here we find simulated systems which do not seem to ever converge to a static state, instead exhibiting repeating cycles of droplet nucleation, elongation of the droplet to a tubular shape by electrostatics and eventually splitting and evaporation of the tube as electrostatic pressure builds up. It is possible that a static state would be reached in a less coarse-grained model applying the electrostatic force separately to the two fluid components as we suggested in discussing Fig. 5b in section III.4, but this remains to be seen.

IV Conclusions

IV.1 Summary

Based on our simulations we developed a theoretical model giving the charge distribution as a function of the fluid composition at each point in time. When antagonicity is of the order of a few $k_{B}T$ this function is linear and for higher antagonicities it is a higher-order polynomial. With this model, we can neglect the complicated dynamics of the ions completely and effectively reduce the system from a quaternary mixture to a unary phase field model. We find that electrostatic interactions in the ion dynamics can in many cases be neglected in the parameter space where mesoscopic structure formation happens, as structure formation occurs only when electrostatic ion fluxes are small in comparison with solvation fluxes. The condition of solvation-dominated ion dynamics $\gamma_{p}>1$, or $\lambda_{I}\gamma<\epsilon\Delta\mu^{2}/6e^{2}$, is identical to the condition of structure formation derived by Onuki Onuki and Kitamura (2004). Assuming for example $\Delta\mu=15k_{B}T$ at a temperature of $T=330$K and $\epsilon_{r}=40$ the product of interface width and surface tension has to fulfill $\lambda_{I}\gamma<10$pN in order for strong fluid-ion coupling and therewith structure formation to be possible. For most water-oil mixtures this value is much larger under normal conditions but can be expected to rapidly decrease to zero as the temperature approaches the critical point of demixing Buhn et al. (2004); Pousaneh et al. (2016). In a recent paper, Okamoto and Onuki have predicted a similar coupling as we observe of charge and order parameter between nonionic solutes in water-oil mixtures to explain the so-called Ouzo-effect Okamoto and Onuki (2018).

By showing equivalency in the case of low antagonicity to the Ohta-Kawasaki model we motivate the observed structure formation and nematic ordering and predict the resulting structure sizes. While it remains to be seen how the scaling laws and morphologies differ from the Ohta-Kawasaki model due to nonlinear coupling of charge and fluid composition for high $\Delta\mu$, we observe excellent agreement with the scaling laws of the Ohta-Kawasaki model up to $\Delta\mu=8k_{B}T.$ Our 3D simulations show essentially the same morphologies as a function of the volume fraction $m$ as is known from the Ohta-Kawasaki model Thomas (2007). On a similar note, Pousaneh and Ciach recently showed in their study of confined mixtures Pousaneh and Ciach (2014) how binary fluid mixtures containing antagonistic salts can also be modeled via the Landau-Brazovskii free energy model, which is recovered by the Ohta-Kawasaki model in the weak-segregation limit Huang and Viñnals (2007).

Our results on nematic ordering point towards the possibility of controlling the average size of ordered domains at a given time or at least the speed of nematic ordering via a number of system parameters affecting $\Lambda$, namely the salt concentration $c_{s}$ but possibly also the temperature, which will strongly affect the ratio of $\lambda_{I}^{3}/\gamma$ close to the critical point. It remains to be studied, whether thermal fluctuations eventually stop nematic ordering or it progresses continuously as suggested by our simulations.

Comparing to experimental data from SANS we numerically reproduce and explain the dissolution of structures at high salt concentration as a result of remixing caused by high electrostatic fluid forces and zero or negative effective surface tension when $\Lambda>\Lambda_{c}.$ For equal volume fractions $m=0.5$ and again assuming $\Delta\mu=15k_{B}T$ and $\epsilon_{r}=40$ this condition is equivalent to $\lambda_{I}^{3}c_{s}^{2}/\gamma>1.9\cdot 10^{26}$N*-1* m*-2* for a critical $\Lambda_{c}\approx 3.1.$ We also interpret the secondary peaks observed in the experimental data as a result of reduced fluid and ion concentrations at the interfaces.

Future work may include quantifying the effects of dielectric permittivity contrasts between the two fluid components, further testing the applicability of the Ohta-Kawasaki model for high, but still realistic values of $\Delta\mu\approx 15k_{B}T,$ a more detailed study of the dynamics of nematic ordering, and reproducing the lamellar phase at low volume fractions, as observed in the experiments of Sadakane et al. Sadakane et al. (2013).

IV.2 Acknowledgements

The authors acknowledge financial support by the Deutsche Forschungsgemeinschaft (DFG) within the Cluster of Excellence "Engineering of Advanced Materials”, the Initiative and Networking Fund of the Helmholtz Association and by the Bavarian Ministry of Economic Affairs and Media, Energy and Technology for the joint projects in the framework of the Helmholtz Institute Erlangen-Nürnberg for Renewable Energy (IEK-11) of Forschungszentrum Jülich. We thank the Jülich Supercomputing Centre and the High Performance Computing Centre Stuttgart for the technical support and the allocated CPU time.