# Field theory description of ion association in re-entrant phase separation
of polyampholytes

Jonas Wessén  Department of Biochemistry, University of Toronto, 1 King’s
College Circle, Toronto, Ontario M5S 1A8, Canada  Tanmoy Pal  Department of
Biochemistry, University of Toronto, 1 King’s College Circle, Toronto, Ontario
M5S 1A8, Canada  Hue Sun Chan  Department of Biochemistry, University of
Toronto, 1 King’s College Circle, Toronto, Ontario M5S 1A8, Canada
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Phase separation of several different overall neutral polyampholyte species
(with zero net charge) is studied in solution with two oppositely charged ion
species that can form ion-pairs through an association reaction. A field
theory description of the system, that treats polyampholyte charge sequence
dependent electrostatic interactions as well as excluded volume effects, is
hereby given. Interestingly, analysis of the model using random phase
approximation and field theoretic simulation consistently show evidence of a
re-entrant polyampholyte phase separation at high ion concentrations when
there is an overall decrease of volume upon ion-association. As an
illustration of the ramifications of our theoretical framework, several
polyampholyte concentration vs ion concentration phase diagrams under constant
temperature conditions are presented to elucidate the dependence of phase
separation behavior on polyampholyte sequence charge pattern as well as ion-
pair dissociation constant, volumetric effects on ion association, solvent
quality, and temperature.

##  I Introduction

Numerous studies in the past decade have indicated that liquid-liquid phase
separation (LLPS) of intrinsically disordered proteins (IDPs) or intrinsically
disordered regions (IDRs) of proteins, often in the presence of nucleic acids
and folded protein domains, is a critical physical mechanism behind the
formation of biologically functional membraneless organelles such as the
nucleolus, cajal bodies and stress granulesBrangwynneHyman2009 ; Rosen12 ;
McKnight12 ; CellBiol ; Nott2015 ; BananiLeeHymanRosen2017 ; ALBERTI2017R1097
; LiChavaliPancsaChavaliBabu ;
MolliexTemirovLeeCoughlinKanagarajKimMittagTaylor ; NatPhys ; Chong2016 ;
Banani2017 ; mingjie2020 . Electrostatics, among other multivalent
interactions, is one of the major drivers of the biological LLPS because
IDP/IDRs generally contain charged amino acid residues in their composition
and are often polyampholytic in nature uversky2002 ; cosb15 ; Robert-Julie ;
SumanPNAS . In addition to a polyampholyte sequence’s charge pattern—or more
generally its charge composition—and concentration SumanJPCB2018 ;
JonasJPCB2021 ; LinPRL ; McCartyDelaneyDanielsenFredricksonShea2019 , the
phase separation propensity of a specific polyampholyte depends on the
condition of the solution determined by such factors as temperature
rolandJACS2019 ; jeetainACS , pressure rolandJACS2019 ; roland20 ;
FetahajJACP2021 , pH WeberJulicher2020 , as well as the concentration and type
of salts present FetahajJACP2021 ; EspinosaJCP2021 . Because of the charged
nature of the polyampholytes, electrostatic interactions from the ion pairs in
the solution affects its LLPS LinPRL ; Kraineretal2020 ; EspinosaJCP2021 ;
FetahajJACP2021  and could be a controlling factor in its bio-engineering
Hong2020 . In general, an oppositely charged pair of ions can stay as solvated
ions or they can form a chemically distinct complex, e.g. solvent-shared ion-
pair, contact ion-pair, through association reactions depending on the
solution condition. Although the effects of ion-strength on polyampholyte LLPS
have been addressed by several analytical and computational studies, the
consequences of ion-association has largely been unexplored.

In general, chemical reactions in biomolecular condensates are of broad
interest because the regulation of biochemical reactions is one of several
major biologically relevant functions of biomolecular
condensatesAlbertiGladfelterMittag2019 . Indeed, several recent computational
studies have addressed phase separation in chemically reactive environments.
For instance, the pH dependence of LLPS was studied by Adame-Arana et al. in a
set-up where the net charge on the phase separating macromolecules is
chemically coupled to the self-ionization of water WeberJulicher2020 . An
investigation by Lin et al. elucidated the role of complex formation between
the SynGAP and PSD-95 molecules in the LLPS of their mixture
LinWuJiaZhangChan2021 . The study by Bartolucci et al. considered both
equilibrium and fuel-driven phase separation of a polymeric component
undergoing an internal molecular transitionBartolucci2021 . These studies
provided valuable insights. However, they treated relevant interactions only
up to mean-field theory (MFT) and did not incorporate amino acid sequence of
the polyampholytes explicitly. Sequence specificity, and generally phase-
separation driven by electrostatic interactions, are inaccessible to the
aforementioned MFT approaches because of the non-neglible contributions from
fluctuationsPopovLeeFredrickson2007 ; LinPRL ; but these effects are
physically and biochemically important. One of the goals of our present work
is to develop theoretical approaches that allow these effects to be tackled.

Pinpointing the exact roles played by all the physical interactions affecting
in vivo LLPS, or even its simplified in vitro counterpart, could be immensely
difficult. For analytical and computational tractability, we consider a simple
model where a polyampholyte species is phase-separating in the presence of two
oppositely charged chemically reacting ion species ${\rm A}^{+}$ and ${\rm
B}^{-}$. We assume that the concentrations of ${\rm A}^{+}$ and ${\rm B}^{-}$
are in thermal equilibrium with the concentration of their charge-neutral
product ${\rm AB}$, following the balance equation

$\ch{A^{+}+B^{-}<=>[K_{\rm a}][K_{\rm d}]AB}.$ (1)

The dissociation/association constants $K_{\rm d}$/$K_{\rm a}$ in Eq. (1) are
defined by

$K_{\rm d}=\frac{1}{K_{\rm a}}=\frac{[{\rm A}^{+}][{\rm B}^{-}]}{[{\rm AB}]},$
(2)

where $[X]$ is the equilibrium concentration of the species $X$ ($={\rm
A}^{+},{\rm B}^{-},{\rm AB}$). The reaction (1) could be used to describe
several chemical processes including self-ionization of water, dissociation of
weak organic acids (e.g. formic acid, acetic acid or carbonic acid) in the
solution, ion association in concentrated solutions of electrolytes at
physical temperatures, at non-polarity solvents or at low temperatures
Valeriani2010 . The dissociation constant $K_{\rm d}$ of a chemical reaction
is an experimentally measurable observable whose value can indicate the
chemical state of the solution. If the initial reactant concentrations are
lower than $K_{\rm d}$, most of the ion pairs are expected to be in the
dissolved state which will result in screening of the polyampholyte’s
electrostatic interactions. On the other hand, when the reactant
concentrations are considerably higher than $K_{\rm d}$, most of the ion pairs
are expected to be in the charge-neutral complex state which will affect the
configuration entropy of the polyampholyte by modulating the effective
excluded volume through steric repulsion. In addition, ion-pair association is
often accompanied by a change in volume SpiroReveszLee1968 ; WarkHsiaSon2008 ;
Marcus1983 . Any such volume change might have an important effect—in addition
to the electrostatic screening effects of the ions—on polyampholyte
conformation at high reactant concentrations. Thus, the chemical state of the
solution determined by Eq. (1) has potential to dictate LLPS behavior of the
polyampholyte species.

At high salt concentrations, non-electrostatic interactions are expected to
play a major role in determining LLPS behavior. Indeed, in a recent explicit
chain molecular dynamics (MD) study, phase separation behaviour of several
proteins at high salt concentrations were attributed to non-electrostatic
interactions such as hydrophobic interactionsKraineretal2020 . In a different
study, the importance of the excluded volume interaction from PEG crowding
agents was highlighted in a system with a phase separating protein that lack
hydrophobic amino acids in its sequence ParkFredrickson2020 . A model that
includes both ion-association along with any volume change upon ion
association and explicit residue level electrostatics of the phase separating
polyampholyte thus offers a unique possibility of capturing many of the
diverse results mentioned above in an unified set-up, at least qualitatively.

With that in mind, here we adopt a trade-off between complexity and
analytical/computational tractability by introducing a simple field theory
model where molecular species in the model interact via excluded volume and
explicit sequence dependent electrostatic interactions. Specifically, we
introduce a bare dissociation constant (corresponding to the dissociation
constant in a solution consisting only of A+, B-, and AB) as a control
parameter for the non-electrostatic energy gain associated with ion-pairing.
To account for the volume change upon association, we introduce further a
relative excluded volume factor $\gamma$ of the product ${\rm AB}$ in Eq. (1)
with respect to the reactants ${\rm A}$ and ${\rm B}$. We study the model
using MFT, random phase approximation (RPA)—where relevant Gaussian level
fluctuations are included above the MFTGonzales-Mozuelos1994 , and fully
fluctuating field theoretic simulation (FTS). We expect the model to be useful
as a base for studying specific systems with suitable modifications.

The structure of this article is as follows. In Sec. II, we introduce our
model and derive its corresponding field theory representation in Sec. II.1.
The model is studied analytically in Sec. III using MFT and RPA, and then
using FTS in Sec. IV. Numerical results obtained from the approximate
analytical calculations and from FTS are shown and compared in Sec. V, and
concluding remarks are given in Sec. VI.

##  II Model definition

Figure 1: (a) Schematic illustration of the effective volume change upon ion-
association, modeled by the parameter $\gamma$. The shading represents
volumetric effect beyond that of the van der Waals volume represented by the
solid outlines of the solutes. (b) Charge sequences considered in this work.
Blue (‘K’) and red (‘E’) beads correspond to charges $+1$ and $-1$,
respectively.

Our system of consideration contains $n_{\rm p}$ linear polymers of identical
composition, each consisting of $N$ residues with electric charges
$\sigma_{\alpha}$, $\alpha=1,\dots,N$. Effects from solvent molecules are
assumed to be implicitly encoded in the microscopic interaction parameters.
All electric charges in this work are given in units of the elementary
protonic charge, and we restrict our calculations to polymers with zero net
electric charge, $\sum_{\alpha=1}^{N}\sigma_{\alpha}=0$. The ‘K’–‘E’ sequences
(Fig. 1 (b)) considered here are representative of the “sv
sequences”DasPappu2013PNAS  used extensively to study LLPS in computational
modelsDasAminLinChan2018 ; PRE2021Pal ; Lin2021 ;
McCartyDelaneyDanielsenFredricksonShea2019 ; SumanJPCB2018 ; JonasJPCB2021 .
An analytically derived single chain property, the sequence charge decoration
(SCD) parameter (defined as
$\mathrm{SCD}=\sum_{\beta=2}^{N}\sum_{\alpha=1}^{\beta-1}\sigma_{\beta}\sigma_{\alpha}\sqrt{\left(\beta-\alpha\right)}/N$),
could be used to discriminate between otherwise charge neutral ‘sv’
sequencesSawleGosh2015 . The phase separation propensity trends of the “sv
sequences” is known to correlate well with their SCD parameter
valuesDasAminLinChan2018 ; McCartyDelaneyDanielsenFredricksonShea2019 ;
JonasJPCB2021 . A ‘K’–‘E’ sequence ($(\mathrm{K}_{8}\mathrm{E}_{8})_{3}$) with
SCD value -4.290 has recently been seen to be phase separating in an in vitro
experimentMadinyaPerry2020 . Compared to that the SCD values of the sequences
studied here are -2.098, -4.349 and -7.374, respectively for sv10, sv15 and
sv20. However, we emphasize that the model presented here is not intended to
be quantitatively accurate given that we ignore some effects that are present
in experimental systems, e.g. the possibility of ion condensation onto the
polymer residues due to attractive interactions beyond the simple Coulomb
forcesSalehi2016 ; Friedowitz2018 ; Ghasemi2020 ; MadinyaPerry2020 ; Sing2020
; Knoerdel2021 . Consequently, while the trend predicted by our theory is
expected to hold for corresponding experimental results, there are
uncertainties NoidReview ; Lin2021  in mapping the theoretical variables such
as temperature to their experimental counterparts for individual systems.

The system further contains equal amounts of unit positive and negative model
ions (denoted ${\rm A}^{+}$ and ${\rm B}^{-}$, respectively) that can undergo
a pair-wise chemical reaction forming a neutral complex ${\rm AB}$ according
to Eq. (1). The number of ${\rm A}$/${\rm B}$ pairs in the dissociated and
bound state, $n_{\rm A}$ ($=n_{\rm B}$) and $n_{\rm AB}$ respectively, are
thus constrained according to

$n_{\rm A}+n_{\rm AB}=n_{\rm T},$ (3)

where $n_{\rm T}$ denotes the total number of ${\rm A}$ and ${\rm B}$ units in
the system.

The position of bead $\alpha$ on polymer $i=1,\dots,n_{\rm p}$ is denoted
${\bf R}_{i,\alpha}$, while the positions of the ${\rm A}^{+}$-, ${\rm
B}^{-}$\- and ${\rm AB}$ particles are represented by $\\{{\bf r}_{{\rm
A},i}\\}_{i=1}^{n_{\rm A}}$, $\\{{\bf r}_{{\rm B},i}\\}_{i=1}^{n_{\rm A}}$ and
$\\{{\bf r}_{{\rm AB},i}\\}_{i=1}^{n_{\rm AB}}$, respectively. Quantities that
depend explicitly on molecule positions are indicated by a hat (e.g.
$\hat{H}$). All energies in this work are given in units of the thermal energy
$k_{\rm B}T$. The total Hamiltonian is
$\hat{H}=\hat{H}_{0}+\hat{H}_{1}+\hat{H}_{2}+\hat{H}_{3}$, where

$\displaystyle\hat{H}_{0}$ $\displaystyle=\frac{3}{2b^{2}}\sum_{i=1}^{n_{\rm
p}}\sum_{\alpha=1}^{N-1}\left({\bf R}_{i,\alpha+1}-{\bf
R}_{i,\alpha}\right)^{2},$ (4) $\displaystyle\hat{H}_{1}$
$\displaystyle=\frac{v}{2}\int\mathrm{d}{\bf r}\hat{\rho}({\bf r})^{2},$
$\displaystyle\hat{H}_{2}$ $\displaystyle=\frac{l_{\rm
B}}{2}\int\mathrm{d}{\bf r}\int\mathrm{d}{\bf r}^{\prime}\frac{\hat{c}({\bf
r})\hat{c}({\bf r}^{\prime})}{|{\bf r}-{\bf r}^{\prime}|}.$

These terms provide, respectively, the energies associated with harmonic chain
connectivity, excluded volume repulsion, and electrostatic interactions. Here
$b$ is the polymer segment length, $v$ is the excluded volume parameter and
$l_{\rm B}=e^{2}/4\pi\epsilon k_{\rm B}T$ is the Bjerrum length ($v$ and
$l_{\rm B}$ modulate the strengths of the excluded volume- and electrostatic
interactions, respectively). We use the segment length $b$ as our unit of
length, and thus values of all dimensionfull quantities (those that carry a
physical unit) are given in powers of $b$. We have expressed $\hat{H}_{1,2}$
in terms of the microscopic particle densities and charge densities, defined
respectively as

$\displaystyle\hat{\rho}({\bf r})$ $\displaystyle=\hat{\rho}_{\rm b}({\bf
r})+\hat{\rho}_{\rm A}({\bf r})+\hat{\rho}_{\rm B}({\bf
r})+\gamma\hat{\rho}_{\rm AB}({\bf r}),$ (5) $\displaystyle\hat{c}({\bf r})$
$\displaystyle=\hat{c}_{\rm b}({\bf r})+\hat{\rho}_{\rm A}({\bf
r})-\hat{\rho}_{\rm B}({\bf r}),$ (6)

where $\hat{\rho}_{\rm b}({\bf r})=\sum_{i,\alpha}\Gamma({\bf r}-{\bf
R}_{i,\alpha})$, $\hat{c}_{\rm b}({\bf
r})=\sum_{i,\alpha}\sigma_{\alpha}\Gamma({\bf r}-{\bf R}_{i,\alpha})$,
$\hat{\rho}_{j}({\bf r})=\sum_{i}\Gamma({\bf r}-{\bf r}_{j,i})$ for $j={\rm
A},{\rm B},{\rm AB}$ and $\gamma$ represents the relative excluded volume
strength of the ${\rm AB}$ complex to a free ion (either ${\rm A}^{+}$ or
${\rm B}^{-}$). Individual polymer beads and ions are treated as Gaussian
distributions $\Gamma({\bf r})=\mathrm{e}^{-{\bf r}^{2}/2a^{2}}/(2\pi
a^{2})^{3/2}$ with smearing length $a=b/\sqrt{6}$ to regulate ultraviolet
divergences arising from contact excluded volume- and electrostatic contact
interactionsWang2010 ; RigglemanRajeevFredrickson2012 . This Gaussian smearing
procedure has also been shown to remedy unphysical binding behaviour in models
of highly concentrated polyelectrolyte solutions that combine equilibrium
constants with Debye-Hückel treatment of electrostatic
interactionsFriedowitz2018 . In this work, we view the Gaussian smearing as a
part of the model definition, rather than as an approximation. While leaving
the long-range behavior of inter-particle forces unaffected, the Gaussian
smearing results in short-range interactions that are considerably “softer”
than typical hard-sphere or Lennard-Jones particle representationsPRE2021Pal .
Consequently, the model is not expected to capture microscopic phenomena such
as strongly oscillating radial distribution functions on short distance scales
that typically characterizes liquid phasesDanielsen2019 .

The parameter $\gamma$ in Eq. (5) may be interpreted through the effective
volume of the ${\rm AB}$ state relative to its dissociated counterpart, with
the sign of $\gamma-2$ corresponding to the sign of ${\rm vol}({\rm AB})-{\rm
vol}({\rm A}^{+})-{\rm vol}({\rm B}^{-})$ as illustrated schematically in Fig.
1 (a). The resulting ${\rm AB}$–${\rm AB}$ and ${\rm AB}$–[other monomer
$={\rm A}^{+}$, ${\rm B}^{-}$, (polymer bead)] excluded volume interaction
strengths are scaled as $\gamma^{2}$ and $\gamma$, respectively, relative to
the [other monomer]–[other monomer] excluded volume repulsion strength $v$. At
the MFT level of an analogous model that takes the ${\rm AB}$ bound state as
an independent system component, these excluded volume interactions result in
effective $\chi$-parameters that are scaled by the same powers of $\gamma$,
and our framework is therefore related to other approaches in the
literatureCruz1995 ; Kudlay_macromolecules_2004 ; Kudlay_JCP_2004  where
effective ion size effects are modeled through appropriate $\chi$-parameter
variations. Implications of effective ion sizes on counter-ion condensation
have also been considered within the context of polyelectrolyte complex
coacervationGhasemi2020 .

The final piece of the Hamiltonian, $\hat{H}_{3}$, models the chemical binding
of the ${\rm AB}$ complex,

$\hat{H}_{3}=(n_{\rm T}-n_{\rm AB})\epsilon_{\rm bound},$ (7)

where $-\epsilon_{\rm bound}$ is the decrement in free energy associated with
a single ${\rm AB}$ binding, in addition to electrostatic and excluded volume
contributions coming from $\hat{H}_{1,2}$. Note that $\epsilon_{\rm bound}$
may include both an energetic component, as well as an entropic component
related e.g. to changes in the number of configuration states, or the entropy
release from reduced electrostriction of nearby solvent moleculesLaidler1987 .
In the remainder of this article, we find it convenient to trade
$\epsilon_{\rm bound}$ in favor of another parameter $K_{0}$, defined as

$K_{0}=\frac{\mathrm{e}^{-\epsilon_{\rm bound}}}{\lambda_{\rm th}^{3}},$ (8)

where $\lambda_{\rm th}$ is the thermal de Broglie wavelength. We show in Sec.
III.1 that $K_{0}$ corresponds to the mean-field dissociation constant in the
dilute limit. This simple physical interpretation of $K_{0}$ breaks down when
particle and charge fluctuations are taken into consideration, as shown in
Sec. III.2.

The partition function of the system,

$Z=\sum_{n_{\rm AB}=0}^{n_{\rm T}}\frac{\int\\{\mathrm{d}{\bf
R}\\}\\{\mathrm{d}{\bf r}_{\rm A}\\}\\{\mathrm{d}{\bf r}_{\rm
B}\\}\\{\mathrm{d}{\bf r}_{\rm
AB}\\}\,\mathrm{e}^{-\hat{H}}}{\tilde{V}^{-(n_{\rm p}+2n_{\rm A}+n_{\rm
AB})}n_{\rm p}!(n_{\rm A}!)^{2}n_{\rm AB}!},$ (9)

is expressed as a sum over the number paired ions $n_{\rm AB}$, where $n_{\rm
A}$ is constrained to $n_{\rm T}-n_{\rm AB}$ according to Eq. (3) and
$\tilde{V}=V/\lambda_{\rm th}^{3}$. In Eq. (9), we use a short-hand notation
$\\{\mathrm{d}{\bf R}\\}\equiv\prod_{i=1}^{n_{\rm
p}}\prod_{\alpha=1}^{N}\mathrm{d}{\bf R}_{i,\alpha}$ and $\\{\mathrm{d}{\bf
r}_{j}\\}\equiv\prod_{i=1}^{n_{j}}\mathrm{d}{\bf r}_{j,i}$ (for $j={\rm
A,B,AB}$) to indicate that the integrals are performed over the positions of
all polymer beads, solvated ions and bound ion pairs. The factorials in the
denominator follow from the fact that the polymers and the ${\rm A}^{+}$,
${\rm B}^{-}$ and ${\rm AB}$ solute are separately indistinguishable. Note
that the particular combination $(n_{\rm A}!)^{2}n_{\rm AB}!$ implies that a
bound state ${\rm AB}$ is distinguishable from any configuration of a free
${\rm A}^{+}$/${\rm B}^{-}$ pair.

The logarithm of the partition function in Eq. (9) can be used to define a
free energy density $f(\rho_{\rm b},\rho_{\rm T};K_{0},l_{\rm
B},v,\dots)=-(\ln Z)/V$, where $V$ is the system volume, $\rho_{\rm b}=n_{\rm
p}N/V$ is the bulk polymer bead number density and $\rho_{\rm T}=n_{\rm T}/V$
is the total bulk number density of A and B ions. We can compute the average
number density of solvated ions $[{\rm A}^{+}]=[{\rm B}^{-}]=\langle n_{\rm
A}\rangle/V$ through the first derivative of $f$ with respect to $K_{0}$,

$[{\rm A}^{+}]=-K_{0}\frac{\partial f}{\partial K_{0}}.$ (10)

The constraint (3) then gives average number density of bound ion pairs as
$[{\rm AB}]=\rho_{\rm T}-[{\rm A}^{+}]$. Knowledge of $[{\rm A}^{+}]$ and
$[{\rm AB}]$ can subsequently be used to compute the dissociation constant
$K_{\rm d}$ according to Eq. (2). Other thermodynamic quantities of interest
to this work are the polymer chain- and ion pair chemical potentials,
$\mu_{\rm p}=N\partial f/\partial\rho_{\rm b}$ and $\mu_{\rm T}=\partial
f/\partial\rho_{\rm T}$ respectively, and the osmotic pressure $\Pi=\rho_{\rm
b}\mu_{\rm p}/N+\rho_{\rm T}\mu_{\rm T}-f$.

Note that $K_{\rm d}$ in Eq. (2) in principle depends on how the bound state
concentration $[{\rm AB}]$ is defined. In this work, we consider the
contributing ${\rm AB}$ states as chemically distint (e.g. through the
appearance of a covalent bond between A and B, an ionic bond between A+ and B-
with strength beyond that of simple electrostatic attraction, or some other
mechanism) from any configuration of free ${\rm A}^{+}$ and ${\rm B}^{-}$
molecules, as follows from the relation in Eq. (10) combined with the
factorials in the denominator in Eq. (9). If $[{\rm AB}]$ is measured by
exploiting the electrostatic properties of the ${\rm A}^{+}$–${\rm B}^{-}$
solution, rather than the chemical nature of the ${\rm AB}$ bond, a different
definition of $[{\rm AB}]$ might be more appropriate where nearby ${\rm
A}^{+}$ and ${\rm B}^{-}$ molecules are counted as effectively bound. This
would be better captured by defining $[{\rm A}^{+}]$ and $[{\rm B}^{-}]$
through their corresponding activity parameters, similarly to how pH is
formally defined through proton activity parameter rather than
concentrationBuck2002pHdef . For simplicity, this alternative approach is not
pursued further in this work.

###  II.1 Deriving the field theory

The basis of the analytical calculations and lattice simulations is a field
representation of the partition function in Eq. (9). To derive its
corresponding field theory, we first decouple the interaction terms in
$\hat{H}_{1}$ and $\hat{H}_{2}$ through standard Hubbard-Stratonovich
transformationsEdwards1965 ; Fredrickson2006 . This introduces two fields
$w({\bf r})$ and $\psi({\bf r})$, conjugate to $\hat{\rho}({\bf r})$ and
$\hat{c}({\bf r})$, respectively, and leads to

$Z=\frac{\int\mathcal{D}w\int\mathcal{D}\psi\,Q_{\rm p}^{n_{\rm
p}}\left(Q_{\rm A}Q_{\rm B}\right)^{n_{\rm
T}}\mathrm{e}^{-H_{0}[w,\psi]}\,S[w,\psi]}{n_{\rm p}!V^{-2n_{\rm
T}}\lambda_{\rm th}^{3n_{\rm T}}K_{0}^{-n_{\rm T}}Z_{w}Z_{\psi}}.$ (11)

up to an inconsequential multiplicative constant, where
$H_{0}[w,\psi]=\int\mathrm{d}{\bf r}[w^{2}/2v+({\bm{\nabla}}\psi)^{2}/8\pi
l_{\rm B}]$. The factor $S[w,\psi]$ contains the sum over $n_{\rm AB}$,

$S=\sum_{n_{\rm AB}=0}^{n_{\rm T}}\frac{1}{n_{\rm AB}![(n_{\rm T}-n_{\rm
AB})!]^{2}}\left(\frac{V^{-1}Q_{\rm AB}}{K_{0}Q_{\rm A}Q_{\rm
B}}\right)^{n_{\rm AB}}.$ (12)

In the above expressions, $Q_{i}\equiv Q_{i}[w({\bf r}),\psi({\bf r})]$ for
$i={\rm p,A,B,AB}$ is the partition function for a type $i$ molecule in
presence of chemical- and electrostatic potential potential fields
$\mathrm{i}\breve{w}({\bf r})\equiv\Gamma\star\mathrm{i}w({\bf r})$ and
$\mathrm{i}\breve{\psi}({\bf r})\equiv\Gamma\star\mathrm{i}\psi({\bf r})$
(where $\star$ denotes spatial convolution, i.e. $\Gamma\star\phi({\bf
r})\equiv\int\mathrm{d}{\bf r}^{\prime}\Gamma({\bf r}-{\bf
r}^{\prime})\phi({\bf r}^{\prime})$ for any function $\phi({\bf r})$). For our
model components, we have

$Q_{\rm p}=\int\mathrm{d}{\bf
R}_{1,\dots,N}\frac{\mathrm{e}^{-\frac{3}{2b^{2}}\sum_{\alpha=1}^{N-1}{\Delta{\bf
R}_{\alpha}}^{2}-\sum_{\alpha=1}^{N}\mathrm{i}\breve{W}_{\alpha}}}{V(2\pi
b^{2}/3)^{3(N-1)/2}},$ (13)

where $\Delta{\bf R}_{\alpha}\equiv{\bf R}_{\alpha+1}-{\bf R}_{\alpha}$,
$\breve{W}_{\alpha}\equiv\breve{w}({\bf
R}_{\alpha})+\sigma_{\alpha}\breve{\psi}({\bf R}_{\alpha})$, and

$\displaystyle Q_{\rm AB}=\frac{1}{V}\int\mathrm{d}{\bf
r}\,\mathrm{e}^{-\gamma\mathrm{i}\breve{w}({\bf r})},\quad Q_{\rm
A}=\frac{1}{V}\int\mathrm{d}{\bf r}\,\mathrm{e}^{-\mathrm{i}[\breve{w}({\bf
r})+\breve{\psi}({\bf r})]},$ (14) $\displaystyle\text{and,}\qquad Q_{\rm
B}=\frac{1}{V}\int\mathrm{d}{\bf r}\,\mathrm{e}^{-\mathrm{i}[\breve{w}({\bf
r})-\breve{\psi}({\bf r})]}.$

Note that all $Q_{i}$ are normalized to $Q_{i}[0,0]=1$. The factors
$Z_{w}\equiv\int\mathcal{D}w\exp(-\int\mathrm{d}{\bf r}w^{2}/2v)$ and
$Z_{\psi}\equiv\int\mathcal{D}\psi\exp(-\int\mathrm{d}{\bf
r}(\nabla\psi)^{2}/8\pi l_{\rm B})$ have been inserted to provide additive
density-independent contributions to the osmotic pressure that cancels the
divergent contributions of the $|{\bf k}|>a^{-1}$ modes of the functional
integrals in Eq. (11)Villet2014 ; Lin2021 .

In the thermodynamic limit, the sum over $n_{\rm AB}$ in $S$ may be replaced
by an integral that can be solved in the saddle-point approximation (this step
is shown in Appendix A). The resulting field picture representation of the
partition function becomes

$Z=\frac{\tilde{V}^{n_{\rm p}+n_{\rm T}}}{n_{\rm p}!n_{\rm
T}!}\frac{\int\mathcal{D}w\int\mathcal{D}\psi\,\mathrm{e}^{-H[w,\psi]}}{Z_{w}Z_{\psi}},$
(15)

where the field Hamiltonian is

$H[w,\psi]=H_{0}-n_{\rm p}\ln Q_{\rm p}-n_{\rm T}\ln Q_{\rm AB}+n_{\rm
T}h_{0}(x).$ (16)

Here, we have introduced a field operator $x[w,\psi]$,

$x[w,\psi]\equiv\frac{2\rho_{\rm T}}{K_{0}}\frac{Q_{\rm AB}}{Q_{\rm A}Q_{\rm
B}}$ (17)

and the function

$h_{0}(x)=\ln\left(1-\frac{\sqrt{1+2x}-1}{x}\right)-\frac{\sqrt{1+2x}-1}{x}.$
(18)

In the following, we find it useful to define two additional functions,
$h_{1}(x)\equiv xh_{0}^{\prime}(x)$ and $h_{2}(x)\equiv xh_{1}^{\prime}(x)$,
where in particular

$h_{1}(x)=\frac{\sqrt{1+2x}-1}{x},$ (19)

satisfies $0\leq h_{1}(x)\leq 1$ for any real $x\geq 0$ (similarly, it may be
shown that $h_{2}(x)\leq 0$). The physical interpretation of the field
operator $h_{1}(x[w,\psi])$ becomes clear when applying the relation in Eq.
(10) to the field partition function in Eq. (15), leading to

$[{\rm A}^{+}]=\langle h_{1}(x)\rangle\rho_{\rm T},$ (20)

i.e. $h_{1}(x[w,\psi])$ is a field operator that averages to the fraction of
dissociated ${\rm A}$/${\rm B}$ pairs. Correspondingly, $[{\rm AB}]=\rho_{\rm
T}\left[1-\langle h_{1}(x)\rangle\right]$ gives the number density of the
bound ion species. The field averaged $h_{1}(x[w,\psi])$ can therefore be used
to evaluate $K_{\rm d}$ in Eq. (2),

$K_{\rm d}=\rho_{\rm T}\frac{\langle h_{1}\rangle^{2}}{1-\langle
h_{1}\rangle}.$ (21)

##  III Analytical theory

To understand the physical implications of our model in Eq. (9), we now
proceed to approximately evaluate the functional integrals in the field
representation of the partition function in Eq. (15). The mean field theory
(MFT) solution, which only accounts for the spatially homogeneous field
configurations, captures the dominant effects on ion dissociation from
excluded volume interactions and, in particular, showcases the crucial role of
the parameter $\gamma$. However, effects from electrostatic interactions are
entirely given by fluctuations of the charge-conjugate field $\psi({\bf r})$,
which are not accounted for in MFTPopovLeeFredrickson2007 . To capture the
leading-order electrostatic effects, we next evaluate Eq. (15) in the random
phase approximation (RPA). This accounts for Gaussian field fluctuations
described by the expansion of $H[w,\psi]$ truncated beyond quadratic order in
$w({\bf r})$ and $\psi({\bf r})$ (corresponding to fluctuations in density and
charge, respectively).

###  III.1 Mean field theory

A spatially homogeneous $\psi({\bf r})=\bar{\psi}$ does not contribute in the
field Hamiltonian in Eq. (16) due to the over-all charge neutrality of the
system, and hence does not contribute in MFT. The MFT solution for $w({\bf
r})=\bar{w}$, given by the vanishing first derivative of $H[w,0]$, satisfies

$\mathrm{i}\bar{w}=v\bar{\rho}_{\rm tot},$ (22)

where

$\bar{\rho}_{\rm tot}=\rho_{\rm
b}+\left[2\bar{h}_{1}+\gamma(1-\bar{h}_{1})\right]\rho_{\rm T}$ (23)

is the total concentration in MFT (counting each AB complex as $\gamma$
units), and $\bar{h}_{1}\equiv h_{1}(x[\bar{w},0])$ with
$x[\bar{w},0]=2\rho_{\rm T}\mathrm{e}^{(2-\gamma)\mathrm{i}\bar{w}}/K_{0}$.
The exponential dependence of $\bar{h}_{1}$ on $\bar{w}$ means that Eq. (22)
generally lacks a closed-form solution, except for in the special case when
$\gamma=2$, and therefore has to be solved numerically for each given set of
parameter values.

Setting $\langle h_{1}\rangle=\bar{h}_{1}$ in Eq. (21) gives the MFT
expression for the dissociation constant, which may be simplified to

$\bar{K}_{\rm d}=K_{0}\exp\left[(\gamma-2)v\bar{\rho}_{\rm tot}\right].$ (24)

In the dilute limit, $\rho_{\rm b},\rho_{\rm T}\rightarrow 0$, this expression
reduces to $\bar{K}_{\rm d}=K_{0}$ verifying the claim in the previous section
that the parameter $K_{0}$ is the dilute-limit MFT dissociation constant. In
the opposite limit, where either $\rho_{\rm b}$ or $\rho_{\rm T}$ are large,
the dissociation constant instead either becomes exponentially suppressed (for
$\gamma<2$) or exponentially enhanced (for $\gamma>2$) at densities
$\gtrsim(|\gamma-2|v)^{-1}$. Physically, this may be interpreted as the bound
${\rm AB}$ state being either favored or disfavored in a dense system
depending on if it yields favorable excluded volume interactions compared to
the dissociated state ${\rm A}^{+}$/${\rm B}^{-}$. At $\gamma=2$, the MFT
dissociation constant is density independent.

The MFT evaluation of the functional integrals
${Z_{w}}^{-1}{Z_{\psi}}^{-1}\int\mathcal{D}w\int\mathcal{D}\psi\,\exp(-H[w,\psi])\approx\exp(-H[\bar{w},0])$
leads to the following expression for the free energy density,

$\bar{f}=-s_{0}+\rho_{\rm T}\bar{h}_{0}+\frac{v}{2}\bar{\rho}_{\rm
tot}\left[\bar{\rho}_{\rm tot}-2(2-\gamma)\bar{h}_{1}\rho_{\rm T}\right],$

where $\bar{h}_{0}\equiv h_{0}(x[\bar{w},0])$, and

$s_{0}=\frac{1}{V}\ln\frac{V^{n_{\rm p}+n_{\rm T}}}{n_{\rm p}!n_{\rm T}!},$

is the configuration entropy density for a system with fully associated ions.
The MFT chemical potentials and osmotic pressure that follow are

$\displaystyle\bar{\mu}_{\rm p}$ $\displaystyle=\ln\frac{\rho_{\rm
b}}{N}+vN\bar{\rho}_{\rm tot},$ $\displaystyle\bar{\mu}_{\rm T}$
$\displaystyle=\ln\left[\rho_{\rm
T}(1-\bar{h}_{1})\right]+v\gamma\bar{\rho}_{\rm tot},$
$\displaystyle\bar{\Pi}$ $\displaystyle=\frac{\rho_{\rm
b}}{N}+(1+\bar{h}_{1})\rho_{\rm T}+\frac{v}{2}\bar{\rho}_{\rm tot}^{2}.$

###  III.2 Gaussian fluctuations:random phase approximation

The Gaussian field fluctuations can be accounted for by expanding the field
Hamiltonian $H[w,\psi]$ in Eq. (16) to quadratic order about the MFT solution,
and then performing the resulting Gaussian functional integrals in the
partition function. This leads to the following correction term to the free
energy density,

$f=\bar{f}+\frac{1}{2}\int\frac{\mathrm{d}{\bf
k}}{(2\pi)^{3}}\ln\left[\frac{4\pi l_{\rm B}v}{{\bf k}^{2}}\det\mathsf{G}({\bf
k})\right],$ (25)

where $\mathsf{G}({\bf k})$ is the 2-by-2 matrix with the Fourier
representation of the coefficients of the quadratic terms in the expansion of
$H[w,\psi]$ and the factor $4\pi l_{\rm B}v/{\bf k}^{2}$ comes from the
product $Z_{w}Z_{\psi}$ in the denominator of Eq. (15). Using the field basis
$(w,\psi)$, we may write

$\mathsf{G}({\bf k})=\mathsf{G}_{0}({\bf k})+\rho_{\rm b}\mathsf{G}_{\rm
p}({\bf k})+\rho_{\rm T}\mathsf{G}_{\rm T}({\bf k}),$

where $\mathsf{G}_{0}({\bf k})={\rm diag}(1/v,{\bf k}^{2}/4\pi l_{\rm B})$,
$\mathsf{G}_{\rm T}({\bf 0})={\rm diag}(-(2-\gamma)^{2}\bar{h}_{2},0)$,
$\mathsf{G}_{\rm T}({\bf k}\neq{\bf 0})=\hat{\Gamma}({\bf k})^{2}{\rm
diag}(\gamma^{2}+(2-\gamma^{2})\bar{h}_{1},2\bar{h}_{1})$, $\mathsf{G}_{\rm
p}({\bf 0})={\rm diag}(0,0)$ and

$\mathsf{G}_{\rm p}({\bf k}\neq{\bf 0})=\hat{\Gamma}({\bf
k})^{2}\begin{pmatrix}g_{\rm mm}({\bf k})&g_{\rm mc}({\bf k})\\\ g_{\rm
mc}({\bf k})&g_{\rm cc}({\bf k})\end{pmatrix}.$

The entries of $\mathsf{G}_{\rm p}({\bf k}\neq{\bf 0})$ are the standard
single chain density-density-, density-charge- and charge-charge correlations
following from the expansion of $Q_{\rm p}$, i.e. $g_{\rm mm}({\bf
k})=\sum_{\alpha,\beta=1}^{N}\mathrm{e}^{-|\alpha-\beta|b^{2}{\bf
k}^{2}/6}/N$, $g_{\rm mc}({\bf
k})=\sum_{\alpha,\beta=1}^{N}\sigma_{\alpha}\mathrm{e}^{-|\alpha-\beta|b^{2}{\bf
k}^{2}/6}/N$ and $g_{\rm cc}({\bf
k})=\sum_{\alpha,\beta=1}^{N}\sigma_{\alpha}\sigma_{\beta}\mathrm{e}^{-|\alpha-\beta|b^{2}{\bf
k}^{2}/6}/N$. In the above expressions, $\hat{\Gamma}({\bf
k})=\mathrm{e}^{-a^{2}{\bf k}^{2}/2}$ is the Fourier transformation of the
Gaussian smearing function $\Gamma({\bf r})$.

The derivatives of Eq. (25) with respect to species numbers, volume and
$K_{0}$ yields corrections from Gaussian field fluctuations to the chemical
potentials, osmotic pressure and fraction of dissociated ion pairs, $\mu_{\rm
p}=\bar{\mu}_{\rm p}+\mu_{\rm p}^{({\rm fl})}$, $\mu_{\rm p}=\bar{\mu}_{\rm
p}+\mu_{\rm p}^{({\rm fl})}$, $\mu_{\rm T}=\bar{\mu}_{\rm T}+\mu_{\rm
T}^{({\rm fl})}$, $\Pi=\bar{\Pi}+\Pi^{({\rm fl})}$ and $\langle
h_{1}(x)\rangle=\bar{h}_{1}+h_{1}^{({\rm fl})}$, respectively. The full
expressions for the RPA corrections $X^{({\rm fl})}$ are given in Appendix B.

The RPA contributions to our thermodynamic observables of interest all involve
integrals over wave numbers ${\bf k}$ that generally need to be computed
numerically. A special case, that can be treated fully analytically, is the
dilute limit of $h_{1}^{({\rm fl})}$, leading to the following RPA expression
for the dilute limit dissociation constant,

$\lim_{\rho_{\rm b},\rho_{\rm T}\rightarrow 0}K_{\rm
d}=K_{0}\left(1+\frac{(2-\gamma^{2})v}{16\pi^{3/2}a^{3}}+\frac{l_{\rm
B}}{\sqrt{\pi}a}\right)^{-1}.$ (26)

This shows that the bare parameter $K_{0}$ can no longer be interpreted as the
dilute limit dissociation constant when fluctuations are included in the
calculation.

##  IV Field Theory Simulations

We complement our findings from RPA calculations with field theoretical
simulations (FTS) that fully capture the fluctuations of the fields $w({\bf
r})$ and $\psi({\bf r})$. Other alternative simulation approaches to
biomolecular LLPS include explicit chain MD simulationsSilmore2017 ;
Dignon2018 ; DasAminLinChan2018 ; SumanJPCB2018  and finite-size scaling
theoryNilssonIrback2020 ; NilssonIrback2021 . In FTS, equilibrium evolution of
the system dictated by the Hamiltonian Eq. (16) is studied by following a
complex Langevin (CL) prescription FredricksonRev2002 ; Fredrickson2006 ;
ParisiWu1981 ; Parisi1983 ; Klauder1983 ; ChanHalpern1986  where the real-
valued continuous fields $w$ and $\psi$ are analytically continued to their
respective complex planes and evolved in the fictitious CL time through the
equations given by

$\displaystyle\frac{\partial w({\bf r},t)}{\partial t}=$
$\displaystyle-\left[\mathrm{i}\tilde{\rho}({\bf r},t)+\frac{w({\bf
r},t)}{v}\right]+\eta_{w}({\bf r},t),$ (27)
$\displaystyle\frac{\partial\psi({\bf r},t)}{\partial t}=$
$\displaystyle-\left[\mathrm{i}\tilde{c}({\bf
r},t)-\frac{\bm{\nabla}^{2}\psi({\bf r},t)}{4\pi l_{\rm
B}}\right]+\eta_{\psi}({\bf r},t).$

This allows ensemble averages to be computed as asymptotic CL time averages.
In Eq. (27), $\eta_{w}$ and $\eta_{\psi}$ are real valued random numbers with
zero mean and variance $2\delta({\bf r}-{\bf
r}^{\prime})\delta(t-t^{\prime})$. The field operators for the number- and
charge densities, $\tilde{\rho}\equiv\tilde{\rho}_{\rm b}+\tilde{\rho}_{\rm
A}+\tilde{\rho}_{\rm B}+\tilde{\rho}_{\rm AB}$ and
$\tilde{c}\equiv\tilde{c}_{\rm b}+\tilde{\rho}_{\rm A}-\tilde{\rho}_{\rm B}$,
respectively, are obtained from

$\tilde{\rho}_{\rm b}({\bf r},t)=\mathrm{i}n_{\rm p}\frac{\delta\ln Q_{\rm
p}}{\delta w({\bf r},t)},\text{ and }\tilde{c}_{\rm b}({\bf
r},t)=\mathrm{i}n_{\rm p}\frac{\delta\ln Q_{\rm p}}{\delta\psi({\bf r},t)},$
(28)

for polymers, and

$\displaystyle\tilde{\rho}_{i}({\bf r},t)$
$\displaystyle=h_{1}(x)\mathrm{i}n_{\rm T}\frac{\delta\ln Q_{i}}{\delta w({\bf
r},t)},\quad i={\rm A,B},$ (29) $\displaystyle\tilde{\rho}_{\rm AB}({\bf
r},t)$ $\displaystyle=\left[1-h_{1}(x)\right]\mathrm{i}n_{\rm
T}\frac{\delta\ln Q_{\rm AB}}{\delta w({\bf r},t)},$

for ions where $x[w,\psi]$ is defined in Eq. (17). Detailed expressions for
the above density operators can be found in Appendix C. Note that
contributions from field fluctuations up to all orders are kept in Eqs. (28)
and (29). We numerically solve Eq. (27) on a cubic lattice with periodic
boundary conditions and lattice spacing $\Delta x$, using a semi-implicit
first order time stepping methodLennon2008SI ; Lin2021  with a CL time-step
$\Delta t=10^{-3}b^{3}$. Use of the semi-implicit time stepping method results
in significantly better numerical stability compared to the Euler–Maruyama
type explicit time-discretization methods Lennon2008SI . As with the CL-time
discretization, the spatial discretization of the continuous fields $w$ and
$\psi$ is an approximation, and, formally, FTS only exactly reproduces the
continuum field theory in the limit $\Delta x\rightarrow 0$. However, the
Gaussian smearing already provides a strong exponential suppression of
contributions from field fluctuations on distance scales smaller than the
smearing length $a$, such that these modes may be omitted with a negligible
numerical effect on physical observables. In this work, we therefore set
$\Delta x=a$ in all FTS computations.

Thermally averaged bulk densities of solvated and bound ion concentrations can
be obtained from the field averaged value of $h_{1}(x)$. Information about how
the components are spatially distributed in the system can be gleaned from
potential of mean forces (PMFs). The PMF $U_{i,j}(r)$ between two components
$i$ and $j$ describe the free-energy landscape for the separation $r$ between
two units of $i$ and $j$, and is related to the corresponding radial
distribution function (RDF) $g_{i,j}(r)$ through the relation

$U_{i,j}(r)=-\ln g_{i,j}(r).$ (30)

Here, $U_{i,j}(r)$ is given in units of the thermal energy $k_{\rm B}T$. Since
explicit particle coordinates have been traded off to the field degrees of
freedom in the field picture, RDFs have to be computed through their field
operators defined by PRE2021Pal

$\displaystyle g_{\rm p,p}(|{\bf r}-{\bf r}^{\prime}|)=$
$\displaystyle\frac{1}{\rho_{\rm b}\rho_{\rm
b}^{\prime}}\bigg{[}\frac{\mathrm{i}}{v}\left\langle\tilde{\rho}_{\rm b}({\bf
r})w({\bf r}^{\prime})\right\rangle-\sum_{i}\left\langle\tilde{\rho}_{\rm
b}({\bf r})\tilde{\rho}_{i}({\bf
r}^{\prime})\right\rangle\bigg{]}-\frac{1}{\rho_{\rm
b}^{\prime}}\frac{\mathrm{e}^{-({\bf r}-{\bf r}^{\prime})^{2}/4a^{2}}}{(4\pi
a^{2})^{3/2}},\quad i=\text{A, B, AB},$ (31) $\displaystyle g_{\rm p,S}(|{\bf
r}-{\bf r}^{\prime}|)=$ $\displaystyle\frac{\left\langle\tilde{\rho}_{\rm
b}({\bf r})\left[\tilde{\rho}_{\rm A}({\bf r}^{\prime})+\tilde{\rho}_{\rm
B}({\bf r}^{\prime})\right]\right\rangle}{2\rho_{\rm b}\rho_{\rm T}\langle
h_{1}\rangle},\quad\text{and},\quad g_{\rm p,B}(|{\bf r}-{\bf
r}^{\prime}|)=\frac{\left\langle\tilde{\rho}_{\rm b}({\bf r})\tilde{\rho}_{\rm
AB}({\bf r}^{\prime})\right\rangle}{\gamma\rho_{\rm b}\rho_{\rm T}(1-\langle
h_{1}\rangle)},$

where the subscripts p, S and B on the RDFs stand for polymer bead, solvated
and bound, respectively. The last term in the expression for $g_{\rm p,p}(r)$
subtracts the contribution from the polymer bead with itself. The RDFs in Eq.
(31) have been normalized to $g_{i,j}=1$ when units of $i$ and $j$ are
uncorrelated, which is expected e.g. at large $r$ if the system contains a
single liquid phase. In Eq. (31), the factor $\rho_{\rm b}^{\prime}=\rho_{\rm
b}-N/V$ provides the correct finite-volume correction to the
polyampholyte–polyampholyte RDF normalisation, and approaches $\rho_{\rm b}$
at large $V$.

If the system is in a globally inhomogeneous state, e.g. by containing several
co-existing macro-phases, certain PMFs may approach non-zero values at large
separations $r$. In particular, the large $r$ behavior of the
polyampholyte–polyampholyte PMF $U_{\rm p,p}(r)$ indicates if the polymers are
homogeneously distributed on large scales or are concentrated in e.g. a dense
droplet. The information of the influences of the solvated ions and the
neutral ion-pairs on phase separation, on the other hand, are obtained from
the polyampholyte–solvated-ions and polyampholyte–ion-pair PMFs, respectively.

##  V Results

Figure 2: The dissociation constant $K_{\rm d}$ (a) and solvated ion density
$[{\rm A}^{+}]=[{\rm B}^{-}]$ (b) computed in RPA (curves) and FTS (circles)
as functions of $\rho_{\rm T}=[{\rm A}^{+}]+[{\rm AB}]$. Black and red
correspond to systems with zero and high polymer densitiy ($\rho_{\rm b}=0$
and $\rho_{\rm b}=5b^{-3}$), respectively. Values of $\gamma$ are indicated by
the linestyles with $\gamma=1.5$, $2$ and $2.5$ corresponding to solid, dashed
and dotted lines, respectively. The gray line in (b) indicates where all ions
exist in the solvated state, i.e. $[{\rm A}^{+}]=\rho_{\rm T}$ and $[{\rm
AB}]=0$. All results are computed using the sv15 charge sequence in Fig. 1(b)
with $l_{\rm B}=0.6b$, $v=0.05b^{3}$, $K_{0}=0.01b^{-3}$ and $a=b/\sqrt{6}$.
The FTS results were computed on $32^{3}$ cubic lattice with lattice spacing
$\Delta x=b/\sqrt{6}$ and a Complex-Langevin time-step $\Delta
t=10^{-3}b^{3}$.

The effects on ion pairing from the solution condition can be understood from
Fig. 2, showing $K_{\rm d}$ and $[{\rm A}^{+}]=[{\rm B}^{-}]$ computed in RPA
and FTS as functions of $\rho_{\rm T}$ at $\gamma=1.5$, $2.0$ and $2.5$ in
polymer empty ($\rho_{\rm b}=0$) and dense ($\rho_{\rm b}=5b^{-3}$) systems
with polyampholyte species sv15. The FTS calculations were done on a $32^{3}$
lattice. The plots in Fig. 2 show excellent agreement between RPA and FTS
except in the very dilute limit where field fluctuations beyond Gaussian order
are expected to become important. Proper sampling of important field
configurations may not be possible through CL methods when the system is
dictated by relatively high excluded volume interaction together with short
range hydrophobicity type interactionsNilsson2022 . The agreement between the
RPA and FTS results seen here is thus reassuring that $\Delta t$ and $\Delta
x$ values used in the FTS implementation provide sufficient numerical
accuracy. When the total density is $\lesssim(v|2-\gamma|)^{-1}$, the dominant
effect on the dissociation constant $K_{\rm d}$ comes from charges in the
surroundings (either free ions or charged polymer residues) that screen out
the electrostatic component of the A-B binding energy, thus increasing the
value of $K_{\rm d}$. In Fig. 2 (a), this effect underlies both the slow
increase of $K_{\rm d}$ with increasing $\rho_{\rm T}$ at $\rho_{\rm b}=0$,
and difference between the $\rho_{\rm b}=0$ and $\rho_{\rm b}=5b^{3}$ curves.
At higher densities, where effects from excluded volume interactions instead
dominate, the exponential dampening (for $\gamma<2$) or enhancement (for
$\gamma>2$) of $K_{\rm d}$ determines whether mainly bound states ${\rm AB}$,
or free ${\rm A}^{+}$ and ${\rm B}^{-}$ are present in the system, as can be
seen in Fig. 2 (b).

Figure 3: RPA phase diagrams in the $(\rho_{\rm b},\rho_{\rm T})$ plane. Sub-
figures (a–c) in the top row correspond to $\gamma=1.5$, $2$ and $2.1$,
respectively, and all with $l_{\rm B}=0.6b$. In the bottom row, sub-figures
(d–f) are for a common $l_{\rm B}=0.9b$ and $\gamma$ values $1.5$, $2$ and
$2.1$, respectively. Other parameter values are $v=0.05b^{3}$,
$K_{0}=0.01b^{-3}$ and $a=b/\sqrt{6}$. The color gradient shows the density of
solvated ions $[{\rm A}^{+}]$, expressed through Eq. (32). Tie-lines indicated
by $\uparrow$ and $\downarrow$ have positive and negative slopes,
respectively.

Fig. 3 shows RPA phase diagrams for the system with sv15 polyampholyte species
computed by matching the osmotic pressure $\Pi(\rho_{\rm b},\rho_{\rm T})$ and
chemical potentials $\mu_{\rm p}(\rho_{\rm b},\rho_{\rm T})$ and $\mu_{\rm
T}(\rho_{\rm b},\rho_{\rm T})$ along the phase boundaryKardar2007 ; Lin2021 .
The binodal curves (solid black curves) enclose the co-existence regions,
where the co-existing bulk density values of $\rho_{\rm b}$ and $\rho_{\rm T}$
are connected by tie-lines (dashed lines). The background color gradient shows
the density of solvated ions, $[{\rm A}^{+}]$, expressed through

${\rm pA}^{+}\equiv-\log_{10}[{\rm A}^{+}]b^{3}\,.$ (32)

The phase diagrams in Fig. 3 display strong dependence on $l_{\rm B}$ and
$\gamma$, and, in particular, predict the possibility of re-entrant phase
separation with increasing $\rho_{\rm T}$. A high $\rho_{\rm T}$ re-entrant
phase separation region occurs at $l_{\rm B}=0.6b$ for $\gamma=1.5$, but
vanishes for $\gamma\geq 2$, which can be understood from the preferred ion-
pairing state in the dense system. As can be seen from Fig. 2, for
$\gamma=1.5$, a large $\rho_{\rm T}$ strongly favors the charge neutral bound
state AB which only interacts through excluded volume repulsion. The resulting
crowding effects leads to an effective attraction between the chains that
further promotes phase separation which enables the re-entrant phase
separation. At $\gamma=2.5$, the high $\rho_{\rm T}$ system contains almost
exclusively solvated ions that strongly screen out the polymer electrostatic
interactions, thus inhibiting phase separation. At the boundary value
$\gamma=2$, although the high $\rho_{\rm T}$ state still contains a
substantial amount of bound ions, the small amount of solvated ions (behaving
roughly as $[A^{+}]\sim\sqrt{\rho_{\rm T}}$) is enough to dissolve the
condensates. When increasing $l_{\rm B}$, the two disconnected co-existence
regions at $\gamma=1.5$, merge into one. Note that the considered
polyampholyte species phase separates for both the $l_{\rm B}$ values used in
Fig. 3 in the absence of any ions.

Figure 4: Slopes of the tielines for the $\gamma=1.5$ and $l_{\rm B}=0.9b$
phase diagram in Fig. 3 (d), as function of the $\rho_{\rm T}$ value of the
polymer-dilute phase. Inset shows the region where tie-line slope changes
sign.

All tie-lines in Fig. 3 are roughly horizontal, but tend to have slightly
positive and negative slopes (indicated by $\uparrow$’s and $\downarrow$’s)
when the ions dominantly are solvated or exist in the bound state,
respectively. This is exemplified in Fig. 4, where we show the tie-line slopes
for the phase diagram in Fig. 3 (d) (i.e. at $l_{\rm B}=0.9b$ for
$\gamma=1.5$). The inset of Fig. 4 shows the region where the tie-line slopes
change sign, beyond which the tie-line slope becomes increasingly more
negative at higher $\rho_{\mathrm{T}}$ values.

Figure 5: Re-entrant phase behaviour dependencies on (a) model parameters and
(b) polyampholyte charge sequece. The black curve in (a) and (b) correspond to
the reference sv15 phase diagram for $\gamma=1.5$ and $l_{\rm B}=0.6b$ in Fig.
3 (top left). Other curves in (a) correspond to varying either $v$, $l_{\rm
B}$ and $K_{0}$. Other curves in (b) instead show the charge sequence
dependence (sequences used are shown in Fig. 1).

In Fig. 5, we focus on the $\gamma=1.5$ and $l_{\rm B}=0.6b$ phase diagram of
Fig. 3 and investigate the dependence of the re-entrant phase behavior on the
model parameters and the polymer charge sequence. The excluded volume
parameter $v$ plays different roles in the upper and lower regions of the re-
entrance phase diagrams, as seen in Fig. 5 (a). While reducing the value of
$l_{\rm B}$ (which reduces the strength of electrostatic interactions) or
increasing $K_{0}$ (giving more solvated ions that provide electrostatic
screening) both decrease the phase separation propensity by shrinking the two
co-existence regions, a reduced $v$ simultaneously enlarges the low $\rho_{\rm
T}$ region and shrinks the high $\rho_{\rm T}$ region. The excluded volume
interactions with ${\rm AB}$ molecules therefore act as to stabilize the high
$\rho_{\rm T}$ condensates, while the excluded volume interactions among the
polymers in the low $\rho_{\rm T}$ region inhibit phase separation. In Fig. 5
(b), we instead swap the sv15 charge sequence by either sv10 or sv20, which
are characterized by smaller or larger blocks of consecutive same-sign
charges, respectively. This degree of “blockiness” can be quantified e.g. by
the $\kappa$ parameter of Das and PappuDasPappu2013 , or by the sequence
charge decoration parameter of Sawle and Gosh SawleGosh2015 , and has been
shown to strongly correlate with phase separation propensityLinPRL ;
SumanJPCB2018 ; DasAminLinChan2018 ;
McCartyDelaneyDanielsenFredricksonShea2019 ; Danielsen2019 . Fig. 5 (b) shows
that the phase diagram exhibits a strong dependence on the charge sequence
blockiness, and that even the topology of the co-existence region may depend
on the sequence (c.f. two disconnected co-existence regions of sv10 and sv15
versus the connected sv20 co-existence region).

Figure 6: FTS results for the species sv15 at $l_{\rm B}=0.6b$, $\rho_{\rm
b}=0.1b^{-3}$, $K_{0}=0.01b^{-3}$and $v=0.05b^{3}$. (a) Representative
snapshots of real non-negative parts of polymer bead density operator
$\tilde{\rho}_{\rm b}$ at various total bulk ion-densities. (b-d) Polymer-
Polymer PMFs at $\gamma=1.5,2$ and $2.5$. (e) and (f) Polymer–Solvated-ions
PMFs and Polymer–ion-pairs PMFs, respetively, at $\gamma=1.5$. The solid lines
and the shaded regions around the solid lines are mean values and standard
deviations of the respective observables computed over five independent
simulations, respectively.

We qualitatively verify some of the LLPS properties obtained from RPA
calculations using FTS calculations of the PMFs $U_{\rm p,p}(r)$, $U_{\rm
p,S}(r)$ and $U_{\rm p,B}(r)$. For the PMF calculations, we now use a larger
($48^{3}$) lattice and focus on the polymer charge sequence sv15. Other fixed
parameter values are $l_{\rm B}=0.6b^{-3}$, $v=0.05b^{3}$, $K_{0}=0.01b^{-3}$
and $\rho_{\rm b}=0.1b^{-3}$. We simulate the system at various initial
reactant concentrations and values of $\gamma$. The results for different
$\rho_{\rm T}$ values are shown in Fig. 6. To obtain reliable statistics, we
ran 5 independent simulations for each $\rho_{\rm T}$. For each independent
run, out of a total $1.5\times 10^{5}$ CL time-steps, we discard the first
$5\times 10^{4}$ steps and use the rest to compute field averages of the RDFs
using a sampling interval of $500$ steps. The equilibration time was taken as
the CL time step after which the real part of the single unit partition
functions stabilized around constant values and the production phase of the
simulations were run until a reasonably small standard deviation of the
computed observables, $K_{\mathrm{d}}$ and PMFs, were reached. In Fig. 6a, we
show representative snapshots of the field picture polymer bead density
$\tilde{\rho}_{\rm b}({\bf r})$, defined in Eq. (28), at different reactant
densities for $\gamma=1.5,2$ and $2.5$. Consistent with the RPA phase diagrams
of Fig. 3, we see polymer droplet formation for all three $\gamma$ values at
relatively low reactant concentrations, but at high reactant concentrations,
droplet formation happens only at $\gamma=1.5$. Different LLPS behavior is
also evident from the polymer-polymer PMF plots shown in Figs. 6b-d where at
the highest reactant concentrations shown, the PMFs go to zero for $\gamma=2$
and $2.5$ and to a small but non-trivial positive value for $\gamma=1.5$ at
large separations. The short-range effective attraction between polymers and
solvated ions or repulsion between polymers and bound ${\rm AB}$ molecules,
respectively, demonstrated by the corresponding PMFs in Fig. 6e-f, are
consistent with the tie-line slopes (see Fig. 4) which tend to be positive and
negative when ions are in their solvated and bound states, respectively.

##  VI Conclusions

In this work, we have introduced a model for polyampholytes undergoing LLPS in
presence of monovalent reactants ${\rm A}^{+}$ and ${\rm B}^{-}$ that can form
a chemically distinct bound state ${\rm AB}$. The dissociated ions have
diametrically different effects on the polymer phase separation from their
electrically neutral bound counterpart: While the free ions screen the LLPS
driving electrostatic interactions between the chains (thus decreasing the
phase separation propensity), the bound pairs ${\rm AB}$ instead function as a
crowding agent that promotes LLPS. Conversely, the crowded and highly charged
environment inside the polyampholyte condensates have non-trivial effects on
the ${\rm A-B}$ dissociation constant compared to the co-existing dilute
phase. This complex interplay between polymer phase separation and ion
dissociation is studied using both analytical and simulation approaches which,
in particular, consistently point towards a novel mechanism for re-entrant
phase behavior under the circumstances where bound ion pairs yield favorable
excluded volume interactions over their dissociated counterpart.

Among the plethora of biological functions of IDP condensates that are
currently being uncovered, regulating chemical reactions seem to be one of
their major functional roles AlbertiGladfelterMittag2019 . Additionally,
chemical reactions in the cellular environment have been shown to be able to
both dissolve condensates and trigger their formation Wippich2013 ;
Zwicker2015 ; Zwicker2017 ; Wurtz2018 ; Kirschbaum2021 . The field theoretic
approach presented in this work constitutes a major methodological advancement
for modelling of such phenomena, due to our treatment of fluctuations compared
to existing MFT theories for polymer phase separation in chemically reactive
environments. Accounting for field fluctuations is necessary to capture amino-
acid sequence dependenceLinPRL , and is essential for describing phase
separation driven by the electrostatic interactionsPopovLeeFredrickson2007
characterized by many IDP species.

The proposed mechanism for ion-triggered re-entrant phase separation relies on
electrostatic screening from solvated ions ${\rm A}^{+}$/${\rm B}^{-}$
combined with crowding effects from their electrically neutral bound state
${\rm AB}$. While Coulombic screening is a long established consequence from
free ions, LLPS promoted by molecular crowding is a relatively less so studied
phenomenonParkFredrickson2020 ; JonasJPCB2021 ; ChanDill1997 ; Ray1971 . Our
model connects these two effects through ion association, quantified by the
bare dissociation constant $K_{0}$. Although it remains to be seen if this
particular mechanism for re-entrance is realized in nature, we nevertheless
believe that our framework for including chemical reactions in RPA and FTS
will be applicable to a wide range of systems through minor phenomenological
modifications.

###### Acknowledgements.

The authors thank Yi-Hsuan Lin and Suman Das for insightful discussions. This
work was supported by Canadian Institutes of Health Research grant NJT-155930
and Natural Sciences and Engineering Research Council of Canada grant
RGPIN-2018-04351 as well as computational resources provided generously by
Compute/Calcul Canada. The data that support the findings of this study are
available from the corresponding authors upon reasonable request. J.W. and
T.P. contributed equally to this work.

##  Appendix A Calculation of the sum over $n_{\rm AB}$ in deriving the field
Hamiltonian

Here, we show how the sum $S[w,\psi]$ over bound ion pairs $n_{\rm AB}$ in Eq.
(12) is computed in the thermodynamic limit. We express $S[w,\psi]$ as

$S=\sum_{n_{\rm AB}=0}^{n_{\rm T}}\frac{y^{n_{\rm AB}}}{n_{\rm AB}![(n_{\rm
T}-n_{\rm AB})!]^{2}},\quad y\equiv\frac{Q_{\rm AB}}{K_{0}VQ_{\rm A}Q_{\rm
B}}.$

In the thermodynamic limit, we can replace the sum by an integral which we
solve in the saddle-point approximation,

$S=\int_{0}^{n_{\rm T}}\mathrm{d}n_{\rm AB}\mathrm{e}^{-I(n_{\rm
AB})}=\mathrm{e}^{-I(n_{\rm AB}^{*})},$

where

$I(n_{\rm AB})=\ln n_{\rm AB}!+2\ln(n_{\rm T}-n_{\rm AB})!-n_{\rm AB}\ln y,$

and $n_{\rm AB}^{*}$ satisfies the saddle condition $I^{\prime}(n_{\rm
AB}^{*})=0$. Using Stirling’s approximation $\ln n!\approx n\ln n-n$ on the
logarithm of the factorials gives the saddle condition

$y=\frac{n_{\rm AB}^{*}}{(n_{\rm T}-n_{\rm AB}^{*})^{2}},$

which has two solutions,

$n_{\rm AB}^{*}=n_{\rm T}\left[1+\frac{1\pm\sqrt{1+2x}}{x}\right],\quad
x=2n_{\rm T}y.$

However, only the ‘$-$’ solution is contained in the integration interval
$0\leq n_{\rm AB}\leq n_{\rm T}$ and contributes to $S[w,\psi]$. Plugging this
solution into $I(n_{\rm AB}^{*})$ gives

$I(n_{\rm AB}^{*})=\ln n_{\rm T}!-n_{\rm T}\ln y+n_{\rm T}h_{0}(x),$

where we again have used Stirling’s approximation in writing $n_{\rm T}\ln
n_{\rm T}-n_{\rm T}\approx\ln n_{\rm T}!$. The function $h_{0}(x)$ is defined
in Eq. (18). The thermodynamic limit of $S[w,\psi]$ then becomes

$S[w,\psi]=\frac{1}{n_{\rm T}!}\left[\frac{Q_{\rm AB}}{K_{0}VQ_{\rm A}Q_{\rm
B}}\right]^{n_{\rm T}}\mathrm{e}^{-n_{\rm T}h_{0}(x)}.$

Plugging this expression into Eq. (11) gives the final field representation of
the partition function in Eq. (15).

##  Appendix B Complete RPA expressions for $\mu_{\rm p,T}$, $\Pi$ and $K_{\rm
d}$

This appendix describes the derivation of the RPA expressions for the chemical
potentials, osmotic pressure and dissociation constant. We start by expanding
the field Hamiltonian in Eq. (16) around the mean field solution $\bar{w}$
given in Eq. (22). To achieve this, we write

$w({\bf r})=\bar{w}+\phi_{0}+\phi({\bf r}),$

with $\int\mathrm{d}{\bf r}\phi({\bf r})\equiv 0$ and $\phi_{0}$ is the ${\bf
k}=\bm{0}$ fluctuation mode. In this field basis, the single molecule
partition functions defined in Eqs. (13) and (14) become

$\displaystyle Q_{\rm p}[\breve{w},\breve{\psi}]=$
$\displaystyle\mathrm{e}^{-\mathrm{i}N(\bar{w}+\phi_{0})}Q_{\rm
p}[\breve{\phi},\breve{\psi}],$ $\displaystyle Q_{\rm
A,B}[\breve{w},\breve{\psi}]=$
$\displaystyle\mathrm{e}^{-\mathrm{i}(\bar{w}+\phi_{0})}Q_{\rm
A,B}[\breve{\phi},\breve{\psi}],$ $\displaystyle Q_{\rm
AB}[\breve{w},\breve{\psi}]=$
$\displaystyle\mathrm{e}^{-\mathrm{i}\gamma(\bar{w}+\phi_{0})}Q_{\rm
AB}[\breve{\phi},\breve{\psi}].$

Charge neutrality results in a global shift symmetry $\psi({\bf
r})\rightarrow\psi({\bf r})+({\rm const.})$ that can be used to eliminate the
${\bf k}=\bm{0}$ mode of $\psi({\bf r})$, such that we can assume
$\int\mathrm{d}{\bf r}\psi({\bf r})=0$ in what follows. The single molecule
partition functions $Q_{i}[\breve{\phi},\breve{\psi}]$ have the quadratic
expansions

$Q_{i}[\breve{\phi},\breve{\psi}]\approx
1-\frac{1}{2V}\int\frac{\mathrm{d}{\bf k}}{(2\pi)^{3}}\Psi^{\rm T}(-{\bf
k})\mathsf{G}_{i}\Psi({\bf k}),$

for $i={\rm p,A,B,AB}$, with $\Psi({\bf k})=\begin{pmatrix}\tilde{\phi}({\bf
k}),\tilde{\psi}({\bf k})\end{pmatrix}^{\rm T}$ containing the Fourier
transformed field fluctuations. Here, $\mathsf{G}_{\rm p}$ is defined as in
the main text and

$\mathsf{G}_{\rm A,B}=\hat{\Gamma}^{2}\begin{pmatrix}1&\pm 1\\\ \pm
1&1\end{pmatrix},\quad\mathsf{G}_{\rm
AB}=\hat{\Gamma}^{2}\begin{pmatrix}\gamma^{2}&0\\\ 0&0\end{pmatrix}.$

To find the RPA expansion of the term $n_{\rm T}h_{0}(x[w,\psi])$ in the field
Hamiltonian, we first express the expansion of the field operator $x[w,\psi]$,
defined in Eq. (17), as

$x[w,\psi]\approx\bar{x}\mathrm{e}^{(2-\gamma)\mathrm{i}\phi_{0}}(1+\epsilon)\equiv
x_{0}(1+\epsilon),$

with $\bar{x}=2\rho_{\rm T}\mathrm{e}^{(2-\gamma)\mathrm{i}\bar{w}}/K_{0}$ and
$\epsilon\equiv\int\mathrm{d}{\bf k}\Psi(-{\bf k})^{\rm T}(\mathsf{G}_{\rm
A}+\mathsf{G}_{\rm B}-\mathsf{G}_{\rm AB})\Psi({\bf k})/2V(2\pi)^{3}$. The RPA
expansion of $h_{0}(x)$ simplifies to $h_{0}(x)\approx
h_{0}(x_{0})+h_{1}(x_{0})\epsilon\approx h_{0}(x_{0})+\bar{h}_{1}\epsilon$
with $\bar{h}_{i}\equiv h_{i}(\bar{x})$. The expansion of $h_{0}(x_{0})$ to
quadratic order in $\phi_{0}$ becomes
$h_{0}(x_{0})\approx\bar{h}_{0}+(2-\gamma)\bar{h}_{1}\mathrm{i}\phi_{0}-(2-\gamma)^{2}\bar{h}_{2}\phi_{0}^{2}$.
This last step can be obtained using the relation
$h_{0}^{\prime\prime}(x)x^{2}=h_{2}(x)-h_{1}(x)$ which can be proven by
considering $(h_{1}(x)x)^{\prime}$ and $h_{i+1}(x)\equiv h^{\prime}_{i}(x)x$.

Combining the above expansion of $h_{0}(x)$ with the expansion of the other
terms in Eq. (16) leads to the RPA expansion of the field Hamiltonian
$H[w,\psi]\approx H[\bar{w},0]+H^{({\rm RPA})}$ with

$H^{({\rm RPA})}=-\frac{n_{\rm
T}}{2}(2-\gamma)^{2}\bar{h}_{2}\phi_{0}^{2}+\frac{1}{2}\int\frac{\mathrm{d}{\bf
k}}{(2\pi)^{3}}\Psi(-{\bf k})^{\rm T}\mathsf{G}\Psi({\bf k})$

where $\mathsf{G}$ is defined as in the main text. One can show that the
contribution from $\phi_{0}$ vanishes in the infinite volume limit
($V\rightarrow\infty$), and will therefore be omitted in the following RPA
calculations. Note, however, that the quadratic expansion of $H[w,\psi]$ is
also used for the semi-implicit CL time integration used in FTS, which
accounts for the ${\bf k}=\bm{0}$ modes explicitly.

Performing the Gaussian field integrals over $\phi({\bf r})$ and $\psi({\bf
r})$ leads to the RPA free energy density in Eq. (25),
$f\approx\bar{f}+f^{({\rm fl})}$ with

$f^{({\rm
fl})}=\frac{1}{4\pi^{2}}\int_{0}^{\infty}\mathrm{d}k\,k^{2}\ln\left(\frac{4\pi
l_{\rm B}v}{k^{2}}\det\mathsf{G}\right).$

The derivatives required for the RPA chemical potentials are

$\frac{\partial f^{({\rm
fl})}}{\partial\rho_{i}}=\frac{1}{4\pi^{2}}\int_{0}^{\infty}\mathrm{d}k\,k^{2}{\rm
Tr}\left[\mathsf{G}^{-1}\frac{\partial\mathsf{G}}{\partial\rho_{i}}\right],\quad
i={\rm b,T}.\\\ $

The density derivatives of the matrix $\mathsf{G}$ includes contributions from
the derivatives of the MFT bound fraction $\bar{h}_{1}$,

$\displaystyle\rho_{\rm T}\frac{\partial\bar{h}_{1}}{\partial\rho_{\rm b}}$
$\displaystyle=\frac{(2-\gamma)\bar{h}_{2}\rho_{\rm
T}v}{1-(2-\gamma)^{2}\bar{h}_{2}\rho_{\rm T}v},$ $\displaystyle\rho_{\rm
T}\frac{\partial\bar{h}_{1}}{\partial\rho_{\rm T}}$
$\displaystyle=\bar{h}_{2}\left[1+(2-\gamma)\rho_{\rm
T}v\frac{\gamma+(2-\gamma)(\bar{h}_{1}+\bar{h}_{2})}{1-(2-\gamma)^{2}\bar{h}_{2}\rho_{\rm
T}v}\right],$

which can be derived by taking the corresponding derivatives of Eq. (22). The
resulting RPA contributions for the chemical potentials, $\mu_{\rm p}^{({\rm
fl})}=N\partial f^{({\rm fl})}/\partial\rho_{\rm b}$ and $\mu_{\rm T}^{({\rm
fl})}=\partial f^{({\rm fl})}/\partial\rho_{\rm T}$, become

$\displaystyle\mu_{\rm p}^{({\rm fl})}$
$\displaystyle=\frac{N}{4\pi^{2}}\int_{0}^{\infty}\mathrm{d}k\,k^{2}{\rm
Tr}\left\\{\mathsf{G}^{-1}\left[\mathsf{G}_{\rm
p}+\frac{(2-\gamma)\bar{h}_{2}\rho_{\rm
T}v}{1-(2-\gamma)^{2}\bar{h}_{2}\rho_{\rm T}v}(\mathsf{G}_{\rm
A}+\mathsf{G}_{\rm B}-\mathsf{G}_{\rm AB})\right]\right\\},$
$\displaystyle\mu_{\rm T}^{({\rm fl})}$
$\displaystyle=\frac{1}{4\pi^{2}}\int_{0}^{\infty}\mathrm{d}k\,k^{2}{\rm
Tr}\left\\{\mathsf{G}^{-1}\left[\mathsf{G}_{\rm
T}+\bar{h}_{2}\left(1+(2-\gamma)\rho_{\rm
T}v\frac{\gamma+(2-\gamma)(\bar{h}_{1}+\bar{h}_{2})}{1-(2-\gamma)^{2}\bar{h}_{2}\rho_{\rm
T}v}\right)(\mathsf{G}_{\rm A}+\mathsf{G}_{\rm B}-\mathsf{G}_{\rm
AB})\right]\right\\}.$

The fluctuation contribution to the osmotic pressure can then be obtained as
$\Pi^{({\rm fl})}=\rho_{\rm b}\mu_{\rm p}^{({\rm fl})}/N+\rho_{\rm T}\mu_{\rm
T}^{({\rm fl})}-f^{({\rm fl})}$.

The RPA expression for the density of dissociated ions $[{\rm A}^{+}]$ follows
from Eq. (10),

$[{\rm A}^{+}]=\bar{h}_{1}\rho_{\rm
T}-\frac{1}{4\pi^{2}}\int_{0}^{\infty}\mathrm{d}k\,k^{2}{\rm
Tr}\left[\mathsf{G}^{-1}K_{0}\frac{\partial\mathsf{G}}{\partial
K_{0}}\right],$

with

$K_{0}\frac{\partial\mathsf{G}}{\partial K_{0}}=-\frac{\rho_{\rm
T}\bar{h}_{2}(\mathsf{G}_{\rm A}+\mathsf{G}_{\rm B}-\mathsf{G}_{\rm
AB})}{1-(2-\gamma)^{2}\rho_{\rm T}v\bar{h}_{2}}.$

The resulting RPA formula for the fraction of dissociated ions becomes
$\langle h_{1}(x)\rangle=\bar{h}_{1}+h_{1}^{({\rm fl})}$ with

$h_{1}^{({\rm fl})}=\frac{\bar{h}_{2}\int_{0}^{\infty}\mathrm{d}k\,k^{2}{\rm
Tr}\left[\mathsf{G}^{-1}(\mathsf{G}_{\rm A}+\mathsf{G}_{\rm B}-\mathsf{G}_{\rm
AB})\right]}{4\pi^{2}(1-(2-\gamma)^{2}\rho_{\rm T}v\bar{h}_{2})},$

which can be plugged into Eq. (21) to obtain the RPA dissociation constant.

##  Appendix C Detailed expressions of the density operators used FTS

For a given field configuration $\\{w({\bf r}),\psi({\bf r})\\}$, the forward
and backward propagators $q_{\mathrm{F}}({\bf r},\alpha)$ and
$q_{\mathrm{B}}({\bf r},\alpha)$, respectively, associated with the discrete
bead-spring polymer chain model used in this work can be used to calculate the
field operators $\tilde{\rho}_{\mathrm{p}}({\bf r})$,
$\tilde{c}_{\mathrm{p}}({\bf r})$ and
$Q_{\mathrm{p}}[\breve{w},\breve{\psi}]$. The chain propagators are
constructed iteratively through Kolmogorov–Chapman equations as

$\displaystyle q_{\mathrm{F}}({\bf r},\alpha+1)$
$\displaystyle=\mathrm{e}^{-\mathrm{i}W({\bf r},\alpha+1)}\Phi\star
q_{\mathrm{F}}({\bf r},\alpha)\,,$ $\displaystyle q_{\mathrm{B}}({\bf
r},\alpha-1)$ $\displaystyle=\mathrm{e}^{-\mathrm{i}W({\bf
r},\alpha-1)}\Phi\star q_{\mathrm{B}}({\bf r},\alpha)\,,$

where $W({\bf r},\alpha)\equiv\breve{w}({\bf
r})+\sigma_{\alpha}\breve{\psi}({\bf r})$, $\Phi({\bf r})\equiv(3/2\pi
b^{2})^{3/2}\exp(-3{\bf r}^{2}/2b^{2})$, and starting from
$q_{\mathrm{F}}({\bf r},1)=\exp[-\mathrm{i}W({\bf r},1)]$ and
$q_{\mathrm{B}}({\bf r},N)=\exp[-\mathrm{i}W({\bf r},N)]$. Here ‘$\star$’
denotes a spatial convolution. One can then show that

$\displaystyle Q_{\mathrm{p}}[\breve{w},\breve{\psi}]$
$\displaystyle=\frac{1}{V}\int\mathrm{d}{\bf r}\,q_{\mathrm{F}}({\bf r},N)\,,$
$\displaystyle\tilde{\rho}_{\mathrm{b}}({\bf r})$
$\displaystyle=\Gamma\star\frac{\rho_{\rm
b}}{N\,Q_{\mathrm{p}}[\breve{w},\breve{\psi}]}\sum_{\alpha=1}^{N}\mathrm{e}^{\mathrm{i}W({\bf
r},\alpha)}q_{\mathrm{F}}({\bf r},\alpha)q_{\mathrm{B}}({\bf r},\alpha)\,,$
$\displaystyle\tilde{c}_{\mathrm{b}}({\bf r})$
$\displaystyle=\Gamma\star\frac{\rho_{\rm
b}}{N\,Q_{\mathrm{p}}[\breve{w},\breve{\psi}]}\sum_{\alpha=1}^{N}\mathrm{e}^{\mathrm{i}W({\bf
r},\alpha)}q_{\mathrm{F}}({\bf r},\alpha)q_{\mathrm{B}}({\bf
r},\alpha)\sigma_{\alpha}\,.$

Bead density operators for $\mathrm{A}^{+}$, $\mathrm{B}^{-}$ and
$\mathrm{AB}$ are given by

$\displaystyle\tilde{\rho}_{\mathrm{A}}$
$\displaystyle=h_{1}(x)\frac{\rho_{\mathrm{T}}}{Q_{\mathrm{A}}[\breve{w},\breve{\psi}]}\mathrm{e}^{-\mathrm{i}\breve{w}-\mathrm{i}\breve{\psi}}\,,$
$\displaystyle\tilde{\rho}_{\mathrm{B}}$
$\displaystyle=h_{1}(x)\frac{\rho_{\mathrm{T}}}{Q_{\mathrm{B}}[\breve{w},\breve{\psi}]}\mathrm{e}^{-\mathrm{i}\breve{w}+\mathrm{i}\breve{\psi}}\,,$
$\displaystyle\tilde{\rho}_{\mathrm{AB}}$
$\displaystyle=\left[1-h_{1}(x)\right]\frac{\gamma\rho_{\mathrm{T}}}{Q_{\mathrm{AB}}[\breve{w},\breve{\psi}]}\mathrm{e}^{-\mathrm{i}\gamma\breve{w}}\,.$

## References

  * [1] Clifford P. Brangwynne, Christian R. Eckmann, David S. Courson, Agata Rybarska, Carsten Hoege, Jöbin Gharakhani, Frank Jülicher, and Anthony A. Hyman.  Germline p granules are liquid droplets that localize by controlled dissolution/condensation.  Science, 324(5935):1729–1732, 2009.
  * [2] P Li, S Banjade, H C Cheng, S Kim, B Chen, L Guo, M Llaguno, J V Hollingsworth, D S King, S F Banani, P S Russ, Q-X Jiang, B T Nixon, and M K Rosen.  Phase transitions in the assembly of multivalent signalling proteins.  Nature, 483:336–340, 2012.
  * [3] M Kato, T W Han, S Xie, K Shi, X Du, L C Wu, H Mirzaei, E J Goldsmith, J Longgood, J Pei, N V Grishin, D E Frantz, J W Schneider, S Chen, L Li, M R Sawaya, D Eisenberg, R Tycko, and S L McKnight.  Cell-free formation of rna granules: low complexity sequence domains form dynamic fibers within hydrogels.  Cell, 149:753–767, 2012.
  * [4] A A Hyman, C A Weber, and F Jülicher.  Liquid-liquid phase separation in biology.  Annu. Rev. Cell Dev. Biol., 30:39–58, 2014.
  * [5] Timothy J Nott, Evangelia Petsalaki, Patrick Farber, Dylan Jervis, Eden Fussner, Anne Plochowietz, Timothy D Craggs, David P Bazett-Jones, Tony Pawson, Julie D Forman-Kay, and Andrew J Baldwin.  Phase transition of a disordered nuage protein generates environmentally responsive membraneless organelles.  Mol. Cell, 57(5):936–947, Mar 2015.
  * [6] Salman F Banani, Hyun O Lee, Anthony A Hyman, and Michael K Rosen.  Biomolecular condensates: organizers of cellular biochemistry.  Nat. Rev. Mol. Cell. Biol., 18:285–298, 02 2017.
  * [7] Simon Alberti.  Phase separation in biology.  Curr. Biol., 27(20):R1097 – R1102, 2017.
  * [8] Xiao-Han Li, Pavithra L. Chavali, Rita Pancsa, Sreenivas Chavali, and M. Madan Babu.  Function and regulation of phase-separated biological condensates.  Biochemistry, 57(17):2452–2461, 2018.
  * [9] Amandine Molliex, Jamshid Temirov, Jihun Lee, Maura Coughlin, Anderson P. Kanagaraj, Hong Joo Kim, Tanja Mittag, and J. Paul Taylor.  Phase separation by low complexity domains promotes stress granule assembly and drives pathological fibrillization.  Cell, 163(1):123 – 133, 2015.
  * [10] C P Brangwynne, P Tompa, and R V Pappu.  Polymer physics of intracellular phase transitions.  Nat. Phys., 11:899–904, 2015.
  * [11] P Andrew Chong and Julie D Forman-Kay.  Liquid–liquid phase separation in cellular signaling systems.  Curr. Opin. Struct. Biol., 41:180–186, 2016.
  * [12] Salman F Banani, Hyun O Lee, Anthony A Hyman, and Michael K Rosen.  Biomolecular condensates: Organizers of cellular biochemistry.  Nat. Rev. Mol. Cell Biol., 18(5):285–298, May 2017.
  * [13] X Chen, X Wu, H Wu, and M Zhang.  Phase transition at the synapse.  Nat. Neurosci., 23:301–310, 2020.
  * [14] V N Uversky.  Natively unfolded proteins: A point where biology waits for physics.  Protein Sci., 11(4):739–756, 2002.
  * [15] T Chen, J Song, and H S Chan.  Theoretical perspectives on nonnative interactions and intrinsic disorder in protein folding and binding.  Curr. Opin. Struct. Biol., 30:32–42, 2015.
  * [16] Robert McCoy Vernon, Paul Andrew Chong, Brian Tsang, Tae Hun Kim, Alaji Bah, Patrick Farber, Hong Lin, and Julie Deborah Forman-Kay.  Pi-pi contacts are an overlooked protein feature relevant to phase separation.  eLife, 7:e31486, feb 2018.
  * [17] Suman Das, Yi-Hsuan Lin, Robert M. Vernon, Julie D. Forman-Kay, and Hue Sun Chan.  Comparative roles of charge, $\pi$, and hydrophobic interactions in sequence-dependent phase separation of intrinsically disordered proteins.  Proceedings of the National Academy of Sciences, 2020.
  * [18] Suman Das, Adam Eisen, Yi-Hsuan Lin, and Hue Sun Chan.  A lattice model of charge-pattern-dependent polyampholyte phase separation.  The Journal of Physical Chemistry B, 122(21):5418–5431, 2018.  PMID: 29397728.
  * [19] Jonas Wessén, Tanmoy Pal, Suman Das, Yi-Hsuan Lin, and Hue Sun Chan.  A simple explicit-solvent model of polyampholyte phase behaviors and its ramifications for dielectric effects in biomolecular condensates.  The Journal of Physical Chemistry B, 125(17):4337–4358, 2021.  PMID: 33890467.
  * [20] Yi-Hsuan Lin, Julie D. Forman-Kay, and Hue Sun Chan.  Sequence-specific polyampholyte phase separation in membraneless organelles.  Phys. Rev. Lett., 117:178101, Oct 2016.
  * [21] James McCarty, Kris T. Delaney, Scott P. O. Danielsen, Glenn H. Fredrickson, and Joan-Emma Shea.  Complete phase diagram for liquid–liquid phase separation of intrinsically disordered proteins.  Mar 2019.
  * [22] S Cinar, H Cinar, H S Chan, and R Winter.  Pressure-sensitive and osmolyte-modulated liquid-liquid phase separation of eye-lens $\gamma$-crystallins.  J. Am. Chem. Soc., 141:7347–7354, 2019.
  * [23] G L Dignon, W Zheng, Y C Kim, and J Mittal.  Temperature-controlled liquid-liquid phase separation of disordered proteins.  ACS Cent. Sci., 5:821–830, 2019.
  * [24] H Cinar, R Oliva, Y-H Lin, X Chen, M Zhang, H S Chan, and R Winter.  Pressure sensitivity of syngap/psd‐95 condensates as a model for postsynaptic densities and its biophysical and neurological ramifications.  Chem Eur. J., 26:11024–11031, 2020.
  * [25] Zamira Fetahaj, Lena Ostermeier, Hasan Cinar, Rosario Oliva, and Roland Winter.  Biomolecular condensates under extreme martian salt conditions.  Journal of the American Chemical Society, 143(13):5247–5259, 2021\.  PMID: 33755443.
  * [26] Omar Adame-Arana, Christoph A. Weber, Vasily Zaburdaev, Jacques Prost, and Frank Jülicher.  Liquid phase separation controlled by ph.  Biophysical Journal, 119(8):1590–1605, Oct 2020.
  * [27] Adiran Garaizar and Jorge R. Espinosa.  Salt dependent phase behavior of intrinsically disordered proteins from a coarse-grained model with explicit water and ions.  The Journal of Chemical Physics, 155(12):125103, 2021.
  * [28] Georg Krainer, Timothy Welsh, Jerelle Joseph, Jorge Espinosa, Ella Csillery, Akshay Sridhar, Zenon Toprakcioglu, Giedre Gudiskyte, Magdalena Czekalska, William Arter, Peter George-Hyslop, Rosana Collepardo-Guevara, Simon Alberti, and Tuomas Knowles.  Reentrant liquid condensate phase of proteins is stabilized by hydrophobic and non-ionic interactions.  05 2020.
  * [29] Kibeom Hong, Daesun Song, and Yongwon Jung.  Behavior control of membrane-less protein liquid condensates with metal ion-induced phase separation.  Nature Communications, 11(1), nov 2020.
  * [30] Simon Alberti, Amy Gladfelter, and Tanja Mittag.  Considerations and challenges in studying liquid-liquid phase separation and biomolecular condensates.  Cell, 176(3):419–434, 2019.
  * [31] Yi-Hsuan Lin, Haowei Wu, Bowen Jia, Mingjie Zhang, and Hue Sun Chan.  Assembly of model postsynaptic densities involves interactions auxiliary to stoichiometric binding.  Biophysical Journal, 2021.
  * [32] Giacomo Bartolucci, Omar Adame-Arana, Xueping Zhao, and Christoph A. Weber.  Controlling composition of coexisting phases via molecular transitions.  Biophysical Journal, 120(21):4682–4697, 2021.
  * [33] Yuri O. Popov, Jonghoon Lee, and Glenn H. Fredrickson.  Field-theoretic simulations of polyelectrolyte complexation.  Journal of Polymer Science Part B: Polymer Physics, 45(24):3223–3230, 2007.
  * [34] Chantal Valeriani, Philip J. Camp, Jos W. Zwanikken, René van Roij, and Marjolein Dijkstra.  Ion association in low-polarity solvents: comparisons between theory, simulation, and experiment.  Soft Matter, 6:2793–2800, 2010.
  * [35] Thomas G. Spiro, Agnes Revesz, and Judith Lee.  Volume changes in ion association reactions. inner- and outer-sphere complexes.  Journal of the American Chemical Society, 90(15):4000–4006, 1968\.
  * [36] Stacey E. Wark, Chih-Hao Hsia, and Dong Hee Son.  Effects of ion solvation and volume change of reaction on the equilibrium and morphology in cation-exchange reaction of nanocrystals.  Journal of the American Chemical Society, 130(29):9550–9555, 2008\.  PMID: 18588299.
  * [37] Y. Marcus.  The volume change on ion-pairing of symmetrical electrolytes.  Zeitschrift für Naturforschung A, 38(2):247–251, 1983.
  * [38] Sohee Park, Ryan Barnes, Yanxian Lin, Byoung-jin Jeon, Saeed Najafi, Kris Delaney, Glenn Fredrickson, Joan Shea, Dong Hwang, and Songi Han.  Dehydration entropy drives liquid-liquid phase separation by molecular crowding.  Communications Chemistry, 3:83, 06 2020.
  * [39] P. González-Mozuelos and M. Olvera de la Cruz.  Random phase approximation for complex charged systems: Application to copolyelectrolytes (polyampholytes).  The Journal of Chemical Physics, 100(1):507–517, 1994.
  * [40] Rahul K. Das and Rohit V. Pappu.  Conformations of intrinsically disordered proteins are influenced by linear sequence distributions of oppositely charged residues.  Proceedings of the National Academy of Sciences, 110(33):13392–13397, 2013.
  * [41] Suman Das, Alan N. Amin, Yi-Hsuan Lin, and Hue Sun Chan.  Coarse-grained residue-based models of disordered protein condensates: utility and limitations of simple charge pattern parameters.  Phys. Chem. Chem. Phys., 20:28558–28574, 2018.
  * [42] Tanmoy Pal, Jonas Wessén, Suman Das, and Hue Sun Chan.  Subcompartmentalization of polyampholyte species in organelle-like condensates is promoted by charge-pattern mismatch and strong excluded-volume interaction.  Phys. Rev. E, 103:042406, Apr 2021.
  * [43] Yi-Hsuan Lin, Jonas Wessén, Tanmoy Pal, Suman Das, and Hue Sun Chan.  Numerical techniques for applications of analytical theories to sequence-dependent phase separations of intrinsically disordered proteins.  Methods in Molecular Biology (Springer-Nature), accepted for publication, arXiv:2201.01920, 2022.
  * [44] Lucas Sawle and Kingshuk Ghosh.  A theoretical method to compute sequence dependent configurational properties in charged polymers and proteins.  The Journal of Chemical Physics, 143(8):085101, 2015.
  * [45] Jason J. Madinya, Li-Wei Chang, Sarah L. Perry, and Charles E. Sing.  Sequence-dependent self-coacervation in high charge-density polyampholytes.  Mol. Syst. Des. Eng., 5:632–644, 2020.
  * [46] Ali Salehi and Ronald G. Larson.  A molecular thermodynamic model of complexation in mixtures of oppositely charged polyelectrolytes with explicit account of charge association/dissociation.  Macromolecules, 49(24):9706–9719, 2016.
  * [47] Sean Friedowitz, Ali Salehi, Ronald G. Larson, and Jian Qin.  Role of electrostatic correlations in polyelectrolyte charge association.  The Journal of Chemical Physics, 149(16):163335, 2018.
  * [48] Mohsen Ghasemi, Sean Friedowitz, and Ronald G. Larson.  Analysis of partitioning of salt through doping of polyelectrolyte complex coacervates.  Macromolecules, 53(16):6928–6945, 2020.
  * [49] Charles E. Sing.  Micro- to macro-phase separation transition in sequence-defined coacervates.  The Journal of Chemical Physics, 152(2):024902, 2020.
  * [50] Ashley R. Knoerdel, Whitney C. Blocher McTigue, and Charles E. Sing.  Transfer matrix model of ph effects in polymeric complex coacervation.  The Journal of Physical Chemistry B, 125(31):8965–8980, 2021.  PMID: 34328340.
  * [51] W. G. Noid.  Perspective: Coarse-grained models for biomolecular systems.  The Journal of Chemical Physics, 139(9):090901, 2013.
  * [52] Zhen-Gang Wang.  Fluctuation in electrolyte solutions: The self energy.  Phys. Rev. E, 81:021501, Feb 2010.
  * [53] Robert A. Riggleman, Rajeev Kumar, and Glenn H. Fredrickson.  Investigation of the interfacial tension of complex coacervates using field-theoretic simulations.  The Journal of Chemical Physics, 136(2):024903, 2012.
  * [54] Scott P. O. Danielsen, James McCarty, Joan-Emma Shea, Kris T. Delaney, and Glenn H. Fredrickson.  Molecular design of self-coacervation phenomena in block polyampholytes.  Proceedings of the National Academy of Sciences, 116(17):8224–8232, 2019.
  * [55] M. Olvera de la Cruz, L. Belloni, M. Delsanti, J. P. Dalbiez, O. Spalla, and M. Drifford.  Precipitation of highly charged polyelectrolyte solutions in the presence of multivalent salts.  The Journal of Chemical Physics, 103(13):5781–5791, 1995.
  * [56] Alexander Kudlay, Alexander V. Ermoshkin, and Monica Olvera de la Cruz.  Complexation of oppositely charged polyelectrolytes: Effect of ion pair formation.  Macromolecules, 37(24):9231–9241, 2004.
  * [57] Alexander Kudlay and Monica Olvera de la Cruz.  Precipitation of oppositely charged polyelectrolytes in salt solutions.  The Journal of Chemical Physics, 120(1):404–412, 2004.
  * [58] Keith J Laidler.  Chemical kinetics.  Number 544.4 LAI. 1987.
  * [59] R. P. Buck, S. Rondinini, A. K. Covington, F. G. K. Baucke, Christopher M. A. Brett, M. F. Camoes, M. J. T. Milton, T. Mussini, R. Naumann, K. W. Pratt, P. Spitzer, and G. S. Wilson.  Measurement of ph. definition, standards, and procedures (iupac recommendations 2002).  Pure and Applied Chemistry, 74(11):2169–2200, 2002.
  * [60] S F Edwards.  The statistical mechanics of polymers with excluded volume.  Proceedings of the Physical Society, 85(4):613–624, apr 1965.
  * [61] Glenn Fredrickson.  The equilibrium theory of inhomogeneous polymers, volume 134.  Oxford University Press on Demand, 2006.
  * [62] Michael C Villet and Glenn H Fredrickson.  Efficient field-theoretic simulation of polymer solutions.  J. Chem. Phys., 141(22):224115, Dec 2014.
  * [63] Kevin S. Silmore, Michael P. Howard, and Athanassios Z. Panagiotopoulos.  Vapour–liquid phase equilibrium and surface tension of fully flexible lennard–jones chains.  Molecular Physics, 115(3):320–327, 2017.
  * [64] Gregory L. Dignon, Wenwei Zheng, Young C. Kim, Robert B. Best, and Jeetain Mittal.  Sequence determinants of protein phase behavior from a coarse-grained model.  PLOS Computational Biology, 14(1):1–23, 01 2018.
  * [65] Daniel Nilsson and Anders Irbäck.  Finite-size scaling analysis of protein droplet formation.  Phys. Rev. E, 101:022413, Feb 2020.
  * [66] Daniel Nilsson and Anders Irbäck.  Finite-size shifts in simulated protein droplet phase diagrams.  The Journal of Chemical Physics, 154(23):235101, 2021.
  * [67] G H Fredrickson, V Ganesan, and F Drolet.  Field-theoretic computer simulation methods for polymers and complex fluids.  Macromolecules, 35:16–39, 2002.
  * [68] G Parisi and Y-S Wu.  Perturbation theory without gauge fixing.  Sci. Sinica, 24:483–496, 1981.
  * [69] G Parisi.  On complex probabilities.  Phys. Lett. B, 131(4):393 – 395, 1983.
  * [70] J R Klauder.  A langevin approach to fermion and quantum spin correlation functions.  J. Phys. A: Math. Gen., 16(10):L317–L319, jul 1983.
  * [71] H S Chan and M B Halpern.  New ghost-free infrared-soft gauges.  Phys. Rev. D, 33:540–547, 1986.
  * [72] Erin M. Lennon, George O. Mohler, Hector D. Ceniceros, Carlos J. García-Cervera, and Glenn H. Fredrickson.  Numerical solutions of the complex langevin equations in polymer field theory.  Multiscale Modeling & Simulation, 6(4):1347–1370, 2008.
  * [73] Daniel Nilsson, Behruz Bozorg, Sandipan Mohanty, Bo Söderberg, and Anders Irbäck.  Limitations of field-theory simulation for exploring phase separation: The role of repulsion in a lattice protein model.  The Journal of Chemical Physics, 156(1):015101, 2022.
  * [74] Mehran Kardar.  Statistical Physics of Particles.  Cambridge University Press, 2007.
  * [75] Rahul K. Das and Rohit V. Pappu.  Conformations of intrinsically disordered proteins are influenced by linear sequence distributions of oppositely charged residues.  Proceedings of the National Academy of Sciences, 110(33):13392–13397, 2013.
  * [76] Frank Wippich, Bernd Bodenmiller, Maria Gustafsson Trajkovska, Stefanie Wanka, Ruedi Aebersold, and Lucas Pelkmans.  Dual specificity kinase dyrk3 couples stress granule condensation/dissolution to mtorc1 signaling.  Cell, 152(4):791–805, 2013.
  * [77] David Zwicker, Anthony A. Hyman, and Frank Jülicher.  Suppression of ostwald ripening in active emulsions.  Phys. Rev. E, 92:012317, Jul 2015.
  * [78] David Zwicker, Rabea Seyboldt, Christoph A Weber, Anthony A Hyman, and Frank Jülicher.  Growth and division of active droplets provides a model for protocells.  Nature Physics, 13(4):408–413, 2017.
  * [79] Jean David Wurtz and Chiu Fan Lee.  Chemical-reaction-controlled phase separated drops: Formation, size selection, and coarsening.  Phys. Rev. Lett., 120:078102, Feb 2018.
  * [80] Jan Kirschbaum and David Zwicker.  Controlling biomolecular condensates via chemical reactions.  Journal of The Royal Society Interface, 18(179):20210255, 2021.
  * [81] Hue Sun Chan and Ken A Dill.  Solvation: how to obtain microscopic energies from partitioning and solvation experiments.  Annual review of biophysics and biomolecular structure, 26(1):425–459, 1997.
  * [82] Ashoka Ray.  Solvophobic interactions and micelle formation in structure forming nonaqueous solvents.  Nature, 231(5301):313–315, 1971.


* [211] Zheng, X., Wise, S., and Cristini, V. Nonlinear simulation of tumor necrosis, neo-vascularization and tissue invasion via an adaptive finite-element/level-set method.  Bulletin of mathematical biology 67, 2 (March 2005), 211—259.