11institutetext: 11email<EMAIL_ADDRESS>22institutetext: 22email:
<EMAIL_ADDRESS>
# Survey on Semantic Stereo Matching / Semantic Depth Estimation

Viny Saajan Victor 11 Peter Neigel 22

###### Abstract

Stereo matching is one of the widely used techniques for inferring depth from
stereo images owing to its robustness and speed. It has become one of the
major topics of research since it finds its applications in autonomous
driving, robotic navigation, 3D reconstruction, and many other fields. Finding
pixel correspondences in non-textured, occluded and reflective areas is the
major challenge in stereo matching. Recent developments have shown that
semantic cues from image segmentation can be used to improve the results of
stereo matching. Many deep neural network architectures have been proposed to
leverage the advantages of semantic segmentation in stereo matching. This
paper aims to give a comparison among the state of art networks both in terms
of accuracy and in terms of speed which are of higher importance in real-time
applications.

###### Keywords:

Stereo Matching, Semantic Segmentation, Depth Estimation, Conditional Random
Fields, Warping Error, Multiscale Context Intertwining

##  1 Introduction

1.1 Stereo Matching: Stereo Matching is a technique used to find point
correspondences between two images of a scene acquired by known cameras. The
traditional approaches such as Semi-Global Matching [5] includes well-defined
steps such as cost matching (calculating matching cost for each pixel in one
of the stereo images depending on the corresponding pixel in the other image),
cost aggregation (calculating matching cost for range of disparity values to
form cost volumes), disparity estimation (calculating the disparity from cost
volume by minimizing the cost function) and disparity post-processing
(improving the disparity). These approaches make assumptions that the scene is
uniformly illuminated and textured.

1.2 Semantic Segmentation: Semantic Segmentation is an image segmentation
methodology that approaches the problem by performing pixel-level
classifications. Traditional approaches include superpixel segmentation [8]
(dividing an image into hundreds of non-overlapping superpixels), active
contour methods [9] (methods which make use of the energy constraints and
forces in the image for separation of region of interest) and watershed
segmentation methods [10] (algorithm which treats pixels values as a local
topography).

1.3 Synergies between Stereo Matching and Semantic Segmentation: Recent works
in semantic matching [6] have shown that the semantic labels obtained from
semantic segmentation improve the accuracy of stereo matching in non-textured,
occluded and reflective regions. Similarly, depth information obtained by
stereo matching can be used to solve possible confusions between similar
semantic categories [3]. This shows that segmentation information and depth
maps complement each other and together they represent high-level information
of the scene.

1.4 State of the Art: Deep learning techniques have shown great success in
both stereo matching and semantic segmentation during the past ten years. The
state of the art uses deep neural networks that are trained for one of these
tasks exclusively and they achieve very reliable accuracy. Recent developments
[6], [3] in the field revealed that stereo matching and semantic segmentation
are intrinsically related, and both pieces of information need to be
considered in an integrated manner to succeed in challenging applications such
as robotics and autonomous navigation. The interactive environment in such
applications requires stereo matching to have a high-level understanding of
the scene. Single deep neural networks jointly trained for both the tasks
achieve better accuracy when compared to standalone models but are not
efficient in terms of speed due to the network complexity. The work presented
in this paper is aimed towards the possible exploration of various methods
used in semantic stereo matching for real-time applications and compare them
in terms of accuracy (network loss) and speed (fps) on benchmark datasets. We
also try to compare, to what extent these methods are utilizing the synergies
between depth information and semantic cues for the benefit of the
applications.

##  2 Semantic Depth Estimation

Semantic Depth Estimation is the method used to improve the depth estimation
obtained from stereo matching using semantic cues(class lael for each pixel in
an image).

###  2.1 Classification

We can majorly classify sematic depth estimation based on two criteria.

Based on training data availability:
Supervised Semantic Depth Estimation: These approaches require datasets with
stereo images and disparity ground truth [1] which are used to train the
network. The network loss is reduced by comparing the output of the network
with the ground truth.

Unsupervised Semantic Depth Estimation: These approaches do not require
disparity ground truth [3], as they mainly rely on warping error. After
calculating the depth and thereby disparity, the right stereo image is warped
with the calculated disparity and compared with the left stereo image. The
network tries to reduce the error as the visual difference between both
images. These approaches still use supervised learning for semantic
segmentation.

Based on network model used:
Standalone Approach: In this approach, networks are trained for depth
estimation and semantic segmentation exclusively and then the result of
semantic segmentation is used to improve disparity.

Joint-Learning Approach: These approaches use a common architecture for both
stereo matching and semantic segmentation in the initial stages to extract
features that are generic for both the tasks and refine the outputs in the
later stages.

Figure 1: Stages involved in semantic depth estimation. Generic features are
extracted from the input images in the joint feature extraction stage. The
output of this stage is fed into 2 networks responsible for disparity
estimation and semantic segmentation. The initial disparity obtained is
improved using semantic cues in disparity refinement stage.

###  2.2 Stages in Extracting Depth from Stereo Images using Semantic Cues

Depth from stereo images using semantic cues can be extracted in 4 generic
stages.

2.2.1 Joint Feature Extraction Stage
This stage extracts generic features from the input stereo images for stereo
matching and semantic segmentation. Features for stereo matching are obtained
in earlier layers and features for semantic segmentation are obtained from
deeper layers since the latter requires more contextual information. The
resolution of the images is usually reduced in this stage, extracting only
relevant features for the task. The output of this stage is given to two
separate subnetworks, one for depth estimation and the other for semantic
segmentation. Hence the network learns a generic representation meaningful for
both the tasks. Dovesi et. al[1] showed that this joint learning is efficient
since a small number of features are sufficient for descent disparity
estimation which in turn significantly increases the frame rate. In [1], the
number of extracted features is given as hyper-parameter which provides
flexibility for the user to choose it according to the requirement (framerate
vs accuracy).

2.2.2 Disparity Estimation Stage
This is the stage where the initial disparity is estimated. The feature map
obtained from the joint feature extraction stage passes through deep
convolution layers to extract features that are specific for disparity
estimation. Cost volumes are created for a range of disparity values and get
the initial disparity map by applying regression. The disparity map obtained
in this stage would contain error in non-textured, occluded and reflective
areas of the images.
Spatial Pyramid Cost Volumes for improving correspondences:

(a)

(b)

Figure 2(a): Construction of spatial pyramid cost volumes from left and right
image features by spatial pooling. Figure is taken from [2]. Figure 2(b):
coarse-to-fine disparity estimation. Each block contains the number of
convolutional layers composing it and the number of features they output and C
is the hyperparameter of the network which corresponds to the number of
extracted features. Figure is taken from [1].

Wu et. al[2] used pyramidal cost volumes to learn the relationship between the
object and its neighbors and thereby acquiring more correspondences. Figure
2(a) shows the construction process of spatial pyramid cost volumes from left
and right image features by spatial pooling. For each resolution of feature
maps, a 4D cost volume is formed by concatenating corresponding unaries from
the left and right image features. The disparity can be inferred from this
volume for that corresponding resolution.
Optimizations to reduce the computation cost keeping desired accuracy: The
network [3] concatenates every feature vector from the left image to all
potential feature vectors from the right image which forms 5-dimensional cost
volumes for both the views. The disparity information from these volumes is
calculated by convolving it with a 3D filter which lets the network learn a
better correlation metric during training. The memory intensity of 3D
convolution is reduced by having an encoder-decoder architecture. In the
network [1], since the feature maps are extracted at different resolutions in
the initial stage (1/16, 1/8, 1/4), this sub-network uses pyramidal network
design where stack of decoders are used to decode the feature map from the
joint feature extractor to estimate coarse-to-fine disparity maps as shown in
Figure 2(b). This stage-wise disparity estimation at different resolution
reduces computation cost as well as allowing the user to manage the speed-
accuracy trade-off dynamically.

2.2.3 Semantic Segmentation Stage
This stage is responsible for extracting semantic labels. Dovesi et. al[1]
used the same coarse-to-fine estimation to achieve the benefits as explained
in the previous stage.
PSP Module [7] for retaining contextual information from different scales: [3]
has a pyramid parsing module for estimating semantic labels. This module
obtains different sub-region representations which are upsampled and
concatenated to form the final feature representation. By doing so both local
and global context information can be carried to the later stage leading to
better accuracy.

2.2.4 Disparity Refinement Stage
Disparity computed in the initial disparity estimation stage contains noise
and its’ accuracy is reduced because of the poor matching in ill-posed
regions. In this stage, disparity results will be refined using semantic cues
based on the assumption that disparity values in the ill-posed regions will be
similar to its neighbors because they belong to the same semantic segment
(local smoothness constraint). The two embeddings (one from disparity
subnetwork and the other from semantic segmentation subnetwork) are
concatenated and this hybrid volume is passed through convolution layers to
get the final refined disparity.
Cost aggregation Module: Multiscale Context Intertwining [11] is used to fuse
the 4D spatial cost volumes from the lowest level to the higher ones in a
recursive way. This module consists of two sub modules.

(a)

(b)

Figure 3: Components of 3D multi-cost aggregation module (a) Hour glass module
is used to upsample the output and (b) 3D feature fusion module is used to
concatenate features. Figure is taken from [2].

i. Hour glass Module: To upsample the output Wu et. al[2] used an hour-glass
model that works on the principle of learnable upsampling. The paper shows
that this approach is efficient because it retains more contextual
information. As shown in figure 3(a), this encoder-decoder architecture
consists of repeated top-down/bottom-up processing. Each module generates a
disparity map for the corresponding resolution.
ii. Feature Fusion Module: Instead of simply concatenating the features, a 3D
feature fusion module is used to retain the context information while fusing
the cost volumes as shown in figure 3(b).
Optimizations to reduce computational cost: During the refinement stage Zhang
et. al[3] focused only on ill-posed regions using residual structures because
mainly those regions are prone to error and thereby reducing computation cost.
[1] compresses the semantic embedding to have the same dimensionality of the
disparity cost volume so that the fusion of these embeddings is computation
effective.

###  2.3 Loss Functions

The overall loss function of the networks [1], [3] and [2] is given by

$\displaystyle L=\alpha L_{disp}+\beta L_{seg}+\gamma L_{ref}$ (1)

where $L_{disp}$ is the initial disparity loss, $L_{seg}$ is the segmentation
loss and $L_{ref}$ is the refined disparity loss,

Loss | Stage | Learning type
---|---|---
[1] L1- Smooth | Initial Disparity Estimation | Supervised
[1] Softmax Cross Entropy | Semantic Segmentation | Supervised
[3] Photometric | Initial Disparity Estimation | Unsupervised
[3] Regularization | Initial Disparity Estimation | Unsupervised
[3] Consistency | Initial Disparity Estimation | Unsupervised
[3] Smoothness | Disparity Refinement | Supervised
[2]Cross-domain discontinuity | Disparity Refinement | Supervised
Table 1: Losses used at different stages in the generic architecture shown in
Figure 1 for Semantic Stereo Matching

$\alpha$, $\beta$ and $\gamma$ are empirically determined constraints.
2.3.1 L1-Smooth Loss [12] : Since disparity estimation from cost volume is a
regression problem [1], [2] use Smooth L1-loss. This loss makes sure that the
estimated disparity is as close as to the ground truth.

$\displaystyle L1_{smooth}=\begin{cases}0.5(d-\hat{d})^{2}&\text{if
$|d-\hat{d}|<1$;}\\\ |d-\hat{d}|-0.5&\text{o}therwise\end{cases}$ (2)

where d is the estimated disparity and $\hat{d}$ is the ground truth
Advantages of using L1-loss are that it has steady gradients for large values
of error and it is less prone to oscillations during updates when the error is
small.
2.3.2 Softmax Cross Entropy [14] : Softmax Cross Entropy loss tries to reduce
the divergence between the predicted probability of the semantic labels and
the ground truth.

$\displaystyle L=-\sum_{c=1}^{M}y_{o,c}\log(p_{o,c})$ (3)

where, $M$ is the total number of classes, $y$ is the binary indicator if
class label $c$ is the correct classification for observation $o$ and $p$ is
the predicted probability observation $o$ is of class $c$.
2.3.3 Photometric Loss [3] : This loss works on the idea that the image
reconstructed from estimated disparity should be very similar to the original
image. It uses Euclidean distance(treats each feature equally) as well as
structure similarity term SSIM $S$ ( metric that quantifies image quality
degradation caused by processing such as data compression or by losses in data
transmission) to improve the robustness in ill-posed regions. The photometric
loss for the left image is

$\displaystyle
L_{p}=\lambda_{1}S(I_{L},I^{{}^{\prime}}_{L})+\lambda_{2}|I_{L}-I^{{}^{\prime}}_{L}|+\lambda_{3}|\nabla
I_{L}-\nabla I^{{}^{\prime}}_{L}|$ (4)

where $I_{L}$ and $I^{{}^{\prime}}_{L}$ are left image and warped right image
respectively, $\nabla$ is the gradient, $\lambda_{1}$,$\lambda_{2}$ and
$\lambda_{3}$ are empirically found constants.
Since the loss depends on the warping error, it doesn’t require a dataset with
ground truth. But this loss introduces high frequency noise.
2.3.4 Regularization Loss [3] : This loss is used to smooth local disparity
from the input images based on the assumption that disparity in the local
region tends to be smooth.

$\displaystyle
L_{r}=\frac{1}{N}\sum|\nabla^{2}_{x}D_{L}|e^{-|\nabla^{2}_{x}I_{L}|}+|\nabla^{2}_{y}D_{L}|e^{-|\nabla^{2}_{y}I_{L}|}$
(5)

where $D_{L}$ is the predicted left disparity, N is the number of pixels,
$\nabla^{2}_{x}$ and $\nabla^{2}_{y}$ are second derivatives along x and y
axes.
This loss reduces the noise introduced from the photometric loss.
2.3.5 Consistency Loss [3] : This loss forces the right and left branches of
the Siamese network to be consistent with each other. The left image can be
warped to the right view and warp it back to the left view by using left
disparity and right disparity. Consistency loss reduces the difference between
this double warped image with the input image.

$\displaystyle
L_{c}=|I_{L}-I^{{}^{\prime\prime}}_{L}|+|I_{R}-I^{{}^{\prime\prime}}_{R}|$ (6)

where $I^{{}^{\prime\prime}}_{L}$ and $I^{{}^{\prime\prime}}_{R}$ are double
warped left and right images respectively.
2.3.6 Smoothness Loss [3] : This loss is used in the refinement stage where
the disparity is improved using semantic cues. This is based on the idea that
disparity should be smooth within a segment. This is conditioned on relatively
good disparity.

$\displaystyle
L_{s}=\sum|\nabla^{2}_{x}D_{L}|(e^{-|\nabla^{2}_{x}f_{L}|}+e^{|Diff-t|})+\nabla^{2}_{y}D_{L}|(e^{-|\nabla^{2}_{y}f_{L}|}+e^{|Diff-t|})$
(7) $\displaystyle Diff=\begin{cases}||D_{L}-D^{{}^{\prime}}_{L}||&\text{if
$||D_{L}-D^{{}^{\prime}}_{L}||\leq t$}\\\ t&\text{o}therwise\end{cases}$ (8)

where $f_{L}$ is the feature vector from the left view and t is an empirically
found threshold.
2.3.7 Cross-domain discontinuity loss [4] : This loss is mainly aimed at
enforcing an explicit link between the two learning tasks by leveraging on the
ground truth pixel-wise semantic labels to improve depth prediction based on
the assumption that depth discontinuities are likely to co-occur with semantic
boundaries.

$\displaystyle
L_{bdry}=\frac{1}{N}\sum_{i,j}(|\varphi_{x}(sem_{i,j})|e^{-|\varphi_{x}(d^{{}^{\prime}}_{i,j})|}+|\varphi_{y}(sem_{i,j})|e^{-|\varphi_{y}(d^{{}^{\prime}}_{i,j})|})$
(9)

where $sem$ is semantic segmentation ground truth, $\varphi_{x}$ and
$\varphi_{y}$ are the intensity gradients between neighboring pixels along the
x and y directions, respectively and $d^{{}^{\prime}}_{i,j}$ is the estimated
disparity between pixels i and j.

##  3 Evaluation

Dataset: KITTI Vision Benchmark Suite 2015 [13] is used as the dataset for the
evaluation which consists of 200 training scenes and 200 test scenes. It also
comprises dynamic scenes for which the ground truth has been established in a
semi-automatic process. It considers a pixel to be correctly estimated if the
disparity or flow end-point error is lesser than 3 pixels.
Metrics used for Evaluation: D1 is the percentage of stereo disparity outliers
in the first frame, all indicate percentage of outliers averaged over all
ground truth pixels, bg indicates the percentage of outliers averaged only
over background regions and fg indicates the percentage of outliers averaged
only over foreground regions.
3.1 Average test set performance: From Table 2, we can infer that [2] gives
the best accuracy on the dataset but it is slower because of the heavy
architecture. [1] gives a real-time performance with comparable accuracy.
Accuracy of [3] is less when compared to the other two networks.

Network | D1-all | Seconds/Frame
---|---|---
RT$S^{2}$ | 3.56 % | 0.02
DispSegNet | 6.33 % | 0.9
SSPCV | 1.85 % | 0.9
Table 2: Shows the average test set performance on KITTI-2015 dataset [13].

$L^{init}_{p}$ | $L^{init}_{c}$ | $L^{init}_{r}$ | $L^{ref}_{p}$ | $L^{ref}_{c}$ | $L^{ref}_{s}$ | $L_{seg}$ | D1-all%
---|---|---|---|---|---|---|---
Yes | Yes | Yes | No | No | No | No | 8.75
Yes | Yes | Yes | Yes | Yes | No | No | 8.60
Yes | Yes | Yes | Yes | Yes | Yes | No | 6.94
Yes | Yes | Yes | Yes | Yes | Yes | Yes | 6.32
Table 3: Ablation study to show the impact of different losses used in [3] on
accuracy.

3.2 Comparison between the methods in each stage
Table 4 shows different modules/methods used in different stages of semantic
depth estimation. From Tables 2 and 4 we can infer that at joint feature
extraction stage, dilated convolutions are giving better accuracy when
compared to encoder-decoder architecture because they preserve the resolution
of the input and spatial pyramid cost volumes are contributing to the increase
in accuracy when compared to single spatial cost volumes beacause they learn
the relationship between object and it’s neighbours to improve
correspondences. Even though hour-glass and 3dff modules used by [2] at the
disparity refinement stage increase the accuracy, they are making the network
slow because of their heavy fusion operations.

Network | JointFeatExt | DispCostVolume | SegCostVolume | DispRefine
---|---|---|---|---
SSPCV | DilatedConvolution | Spatial Pyramid | Single Spatial | Hour-Glass and 3DFF
RTS2 | Encoder architecture | Spatial Pyramid | Spatial Pyramid | CascadeResidualConcatenate
DispSegNet | ResNet50(Siamese) | Single Spatial | Single Spatial | Concatenation
Table 4: Different methods/modules used in semantic depth estimation stages

3.3 Ablation Studies:
In this section, we compare how much each module in the network is
contributing to the overall objective.
DispSegNet: Table 3 shows the ablation study for different loss terms in [3]
where, $init$ is the initial disparity estimation, $ref$ is the disparity
refinement, p is the photometric loss, c is consistency loss, r is
regularization loss, s is smoothness loss and $seg$ is the supervised semantic
segmentation loss. From the Table, we can observe that smoothness loss in the
refinement stage is contributing significantly to the accuracy because they
refine the disparity using semantic cues.
SSPCV: We can observe from Table 5 that pyramid cost volumes, feature fusion
module are contributing to the accuracy beacause of context information
retention and boundary loss function is predominantly contributing to the
accuracy since it is aimed at enforcing an explicit link between the two
learning tasks by leveraging on the ground truth pixel-wise semantic labels to
improve depth prediction.
RTS2: We observe from Table 6 that stage 2 of [1] (i.e 1/8 resolution) is
giving comparable accuracy with very good framerate.

SemBranch | PyramidCost | DilatedConv | FFM | BoundaryLoss | D1-all%
---|---|---|---|---|---
Yes | No | No | No | No | 2.37
Yes | Yes(3DMultiCost) | No | No | No | 1.99
Yes | Yes | Yes | No | No | 2.10
Yes | Yes | Yes | Yes | No | 1.93
Yes | Yes | Yes | Yes | Yes | 1.85
Table 5: Ablation study to compare different model variants of [2] to justify the selection criterion. Init Disp Module | Sem Seg Module | Disp Refine Module | Stage | D1-all% | FPS
---|---|---|---|---|---
Yes | Yes | Yes | 1 | 8.00 | 17.2
Yes | Yes | Yes | 2 | 4.70 | 10.9
Yes | No | No | 3 | 3.98 | 8.1
Yes | Yes | No | 3 | 3.90 | 6.6
Yes | Yes | Yes | 3 | 3.33 | 6.3
Table 6: Ablation study to show the impact of different modules of [1] on
accuracy and framerate.

##  4 Conclusion

From all the comparisons between the networks, we can infer that the networks
which are focusing on accuracy are slower because of their heavy architecture.
And networks focusing on framerate are compromising on accuracy. One has to
always decide on the trade-off between time and accuracy depending on the
application requirements. We also examined network which provides the
flexibility for the user to choose between these two entities depending on
their requirements. As future improvement, instance segmentation can be used
for improving the depth estimation.

## References

  * [1] Dovesi, Pier Luigi and Poggi, Matteo and Andraghetti, Lorenzo and Martí, Miquel and Kjellström, Hedvig and Pieropan, Alessandro and Mattoccia, Stefano, “Real-time semantic stereo matching, 2020 IEEE International Conference on Robotics and Automation (ICRA), IEEE, pp. 10780–10787, 2020.
  * [2] Wu, Zhenyao and Wu, Xinyi and Zhang, Xiaoping and Wang, Song and Ju, Lili, “Semantic stereo matching with pyramid cost volumes, Proceedings of the IEEE/CVF International Conference on Computer Vision, IEEE, pp. 7484–7493, 2019.
  * [3] Zhang, Junming and Skinner, Katherine A and Vasudevan, Ram and Johnson-Roberson, Matthew, “Dispsegnet: Leveraging semantics for end-to-end learning of disparity estimation from stereo imagery, IEEE Robotics and Automation Letters, IEEE, vol. 4, no. 2, pp. 1162–1169, 2019.
  * [4] Ramirez, Pierluigi Zama and Poggi, Matteo and Tosi, Fabio and Mattoccia, Stefano and Di Stefano, Luigi, “Geometry meets semantics for semi-supervised monocular depth estimation, Asian Conference on Computer Vision, Springer, pp. 298–313, 2018.
  * [5] Hirschmuller, Heiko, “Accurate and efficient stereo processing by semi-global matching and mutual information, 2005 IEEE Computer Society Conference on Computer Vision and Pattern Recognition (CVPR’05), IEEE, vol. 2, pp. 807–814, 2005.
  * [6] Yang, Guorun and Zhao, Hengshuang and Shi, Jianping and Deng, Zhidong and Jia, Jiaya, “Segstereo: Exploiting semantic information for disparity estimation, Proceedings of the European Conference on Computer Vision (ECCV), pp. 636–651, 2018.
  * [7] Zhao, Hengshuang and Shi, Jianping and Qi, Xiaojuan and Wang, Xiaogang and Jia, Jiaya, ‘Pyramid scene parsing network, Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 2881–2890, 2017.
  * [8] Zhang, Yongxia and Li, Xuemei and Gao, Xifeng and Zhang, Caiming, “A simple algorithm of superpixel segmentation with boundary constraint, IEEE transactions on circuits and systems for video technology, IEEE, vol. 27, no. 7, pp. 1502–1514, 2016.
  * [9] Xie, Xianghua, “Textured image segmentation using active contours, International Conference on Computer Vision, Imaging and Computer Graphics, Springer, pp. 357–369, 2009.
  * [10] Preim, Bernhard and Botha, Charl, “Image analysis for medical visualization. Visual computing for medicine, Elsevier, 2014.
  * [11] Lin, Di and Ji, Yuanfeng and Lischinski, Dani and Cohen-Or, Daniel and Huang, Hui, “Multi-scale context intertwining for semantic segmentation, Proceedings of the European Conference on Computer Vision (ECCV), pp. 603–619, 2018.
  * [12] Huber, Peter J, “Robust estimation of a location parameter, Breakthroughs in statistics, Springer, pp. 492–518, 1992.
  * [13] Menze, Moritz and Geiger, Andreas, “Object scene flow for autonomous vehicles, Proceedings of the IEEE conference on computer vision and pattern recognition, pp. 3061–3070, 2015.
  * [14] Goodfellow, Ian and Bengio, Yoshua and Courville, Aaron, “Deep Learning, IEEE Robotics and Automation Letters, MIT press, 2016.


# Exact Analytical Solution of the Flory-Huggins Model and Extensions to
Multicomponent Systems

J. Pedro de Souza  Omenn-Darling Bioengineering Institute, Princeton
University, Princeton, NJ, 08544  Howard A. Stone  Department of Mechanical
and Aerospace Engineering, Princeton University, Princeton, NJ, 08544

###### Abstract

The Flory-Huggins theory describes the phase separation of solutions
containing polymers. Although it finds widespread application from polymer
physics to materials science to biology, the concentrations that coexist in
separate phases at equilibrium have not been determined analytically, and
numerical techniques are required that restrict the theory’s ease of
application. In this work, we derive an implicit analytical solution to the
Flory-Huggins theory of one polymer in a solvent by applying a procedure that
we call the implicit substitution method. While the solutions are implicit and
in the form of composite variables, they can be mapped explicitly to a phase
diagram in composition space. We apply the same formalism to multicomponent
polymeric systems, where we find analytical solutions for polydisperse
mixtures of polymers of one type. Finally, while complete analytical solutions
are not possible for arbitrary mixtures, we propose computationally efficient
strategies to map out coexistence curves for systems with many components of
different polymer types.

##  I Introduction

Since its introduction in the early 1940s first by Huggins Huggins (1941) and
then by Flory Flory (1942), the Flory-Huggins theory has become the most
widespread model of phase separation of polymer solutions. Due to its
simplicity in construction, the theory, with a minimal set of parameters,
captures the main trends in the polymer partitioning between phases. The
theory has been applied in numerous contexts including in chemical processing
Cantow (1967); Barton (2018), materials science Milczewska and Voelkel (2006);
Peinemann  _et al._ (2007); Angelopoulou  _et al._ (2023); Callaway  _et al._
(2018), and also in the burgeoning field of biomolecular condensation
Brangwynne  _et al._ (2015); Berry  _et al._ (2018); Pappu  _et al._ (2023).
While many corrections and modifications have been proposed to capture more
detailed physical phenomena Overbeek and Voorn (1957); Simha and Somcynsky
(1969); Scheutjens and Fleer (1980); Tanaka (1989); Tanaka and Stockmayer
(1994); Semenov and Rubinstein (1998); Ruzette and Mayes (2001); Wang (2008,
2017); Wessén  _et al._ (2021), the original Flory-Huggins theory is the
standard starting point for any analysis of materials that include polymeric
components.

From its physics-based construction, the Flory-Huggins theory outputs a system
of nonlinear constraints specifying thermodynamic equilibria that define the
coexisting phase concentrations. While the equilibrium conditions can be
solved numerically, up to this point, to our knowledge, no exact analytical
solution has been found. Therefore, any evaluation of the theory or fitting of
the theory to data requires numerical solutions, hampering the ease of use of
the theory. In some regimes, analytical approximations can be applied Sanchez
(1984); Tsenoglou and Papaspyrides (2001); Scheinhardt-Engels  _et al._
(2004); Qian  _et al._ (2022); van Leuken  _et al._ (2023a), but these
approximations may involve complicated expressions and may not be uniformly
valid, even for a single polymer-solvent combination.

The technical challenge of describing coexisting phases is more apparent for
multicomponent mixtures Mao  _et al._ (2019); Jacobs and Frenkel (2017);
Jacobs (2023), even for mixtures of the same polymer type with different
polymer lengths Huggins and Okamoto (1967). In order to construct a phase
diagram using direct numerical solution, it is necessary to discretize the
compositional space, and solve nonlinear systems of equations at each point.
As the number of polymer types, $M$, increases, the discretized space
increases exponentially, with a power proportional to $M$, which corresponds
to a large number of discrete compositions, especially if fine resolution is
desired. In many applications, such as in cellular biology, the number of
polymer types can exceed tens of thousands Wühr  _et al._ (2014). Such a brute
force numerical approach may not be feasible to describe the breadth of
coexisting concentrations over the full compositional space for $M>4$, never
mind for $M>10^{4}$.

Figure 1: Schematic of the systems under study, their corresponding variables,
and the solution strategy employed in each case for (a) a single polymer-
solvent solution, (b) a polydisperse polymer-solvent solution, and (c) a
solution containing multiple polymer types in solvent.

Here, we analytically derive an exact implicit solution to the Flory-Huggins
theory of a one polymer-one solvent system using an implicit substitution
method. The key idea is to solve for each constraint for the Flory-Huggins
$\chi$ parameter, and then substitute the value of $\chi$ into each equation
so as to eliminate it from the system. We explore the solution’s asymptotic
behavior, referring back to the original applications of the theory by Flory
Flory (1942). The general idea of our approach for this and other problems
described in this paper is illustrated in Figure 1.

Next, we apply the implicit substitution method to multicomponent mixtures
where the relative magnitude of the pair interactions are known. There, we
discover an exact analytical solution for polydisperse polymer samples of the
same type, and distill the phase diagram of these mixtures into two
independent variables that fully specify the coexistence curves. Using this
strategy, we can analytically describe the composition (molecular weight
distribution and total polymer volume fraction) of the coexisting phases.

For generalized mixtures, we use the implicit substitution method to simplify
a nonlinear system of $M+1$ equations down to one nonlinear “master” equation
with one unknown— the solvent partitioning between the phases. While this
reduction leads to great simplification of the equilibrium calculation, in the
case of arbitrary mixtures, the master equation must be solved numerically. If
desired, we can discretize the compositional space and construct a $\chi$
surface by solving just this one equation at each point—and we demonstrate
this idea for a two polymer, one solvent mixture. To minimize the number of
function evaluations, we propose a sampling method for the implicit function
that obviates the need to discretize the entire compositional space. Choosing
sampling points only requires finding the poles of the implicit $\chi$
function at a fixed distance from the global critical point, which we
demonstrate, although we leave a full computational implementation of the
compositional sampling for a large number of components to future work.

The implicit substitution method developed in this work is a powerful
technique to solve or simplify the nonlinear systems of equations encountered
in the thermodynamics of solutions. The only price to be paid is that the
derived solutions are implicit in composite composition variables, themselves
a function that is a combination of the coexisting phase concentrations.
Nevertheless, these composite composition variables are experimentally
accessible, and we expect the technique will find applications in
thermodynamic analyses in various contexts.

##  II The Flory-Huggins Model: One Polymer

###  II.1 Theoretical background

The Flory-Huggins model may be derived using a free energy density where a
polymer composed of $N$ monomers occupies a lattice. The free energy density
per lattice site, $f$, may be expressed in dimensionless terms, $\tilde{f}$,
as:

$\begin{split}\tilde{f}=\frac{fv}{k_{B}T}=&\frac{1}{N}\phi\ln(\phi)+(1-\phi)\ln(1-\phi)\\\
&+\chi\phi(1-\phi).\end{split}$ (1)

Here, $v$ is a lattice site volume, $k_{B}$ is the Boltzmann constant, $T$ is
the absolute temperature, and $\phi$ is the polymer volume fraction. $\chi$ is
the so-called Flory parameter, which captures the interaction between solvent
and the polymer Doi (2013).

From this free energy, we can define the dimensionless chemical potential of
the polymer, ${\mu}=\partial\tilde{f}/\partial\phi$, up to an arbitrary
constant,

${\mu}=\frac{1}{N}\ln(\phi)-\ln(1-\phi)-2\chi\phi,$ (2)

and the dimensionless osmotic pressure of the solution is
$\tilde{\Pi}=-\tilde{f}+\phi\tilde{f}^{\prime}$, so that

${\Pi}=\left(\frac{1}{N}-1\right)\phi-\ln(1-\phi)-\chi\phi^{2}.$ (3)

The model predicts two coexisting phases, a dense phase and a dilute phase of
the polymer, above a critical $\chi$ value, $\chi>\chi_{c}$,

$\chi_{c}=\frac{1}{2}\left(1+\frac{1}{\sqrt{N}}\right)^{2}$ (4)

and beginning at a critical polymer volume fraction $\phi=\phi_{c}$:

$\phi_{c}=\frac{1}{1+\sqrt{N}}.$ (5)

It is found that as $N$ increases, the phase diagram becomes asymmetric
relative to the center line at $\phi=1/2$.

Chemical equilibrium between two phases $A$ and $B$ at volume fractions
$\phi_{A}$ and $\phi_{B}$ is specified by two contraints—equal chemical
potential of the polymer and equal osmotic pressure of the solvent in the two
phases:

$\mu(\phi_{A})=\mu(\phi_{B})$ (6a) $\Pi(\phi_{A})=\Pi(\phi_{B}).$ (6b)

These coupled equations can be rewritten explicitly as:

$\frac{1}{N}\ln\left(\frac{\phi_{A}}{\phi_{B}}\right)-\ln\left(\frac{1-\phi_{A}}{1-\phi_{B}}\right)-2\chi(\phi_{A}-\phi_{B})=0$
(7a)
$\left(\frac{1}{N}-1\right)(\phi_{A}-\phi_{B})-\ln\left(\frac{1-\phi_{A}}{1-\phi_{B}}\right)-\chi(\phi_{A}^{2}-\phi_{B}^{2})=0.$
(7b)

If $\chi$ is known, then we have two equations with two unknowns, $\phi_{A}$
and $\phi_{B}$. To our knowledge, these coupled equations are impossible to
deconvolve into two explicit functions describing $\phi_{A}$ and $\phi_{B}$
with known analytical techniques. However, in order to make analytical
progress, we propose two simultaneous strategies. First, we search for
implicit solutions, where we treat $\chi$ as an unknown, so as to gain an
additional degree of freedom in the compositional space; also, we can
eliminate $\chi$ from the system of equations quite easily. Second, instead of
working with the variables $\phi_{A}$ and $\phi_{B}$, we use two orthogonal
composite variables that combine $\phi_{A}$ and $\phi_{B}$, which we call $y$
and $z$, that become apparent after eliminating $\chi$. Thus, the problem is
reduced to solving for one unknown variable, $z$, in terms of one known
variable $y$. We call this solution method an implicit substitution method,
since we use substitution to eliminate the implicit function $\chi$. Now that
the strategy has been explained, we implement it in detail below.

###  II.2 Analytical solution

First, we solve the equations (7a) and (7b) for $\chi$:

$\chi=\frac{\frac{1}{N}\ln\left(\frac{\phi_{A}}{\phi_{B}}\right)-\ln\left(\frac{1-\phi_{A}}{1-\phi_{B}}\right)}{2(\phi_{A}-\phi_{B})}$
(8a)
$\chi=\frac{\left(\frac{1}{N}-1\right)(\phi_{A}-\phi_{B})-\ln\left(\frac{1-\phi_{A}}{1-\phi_{B}}\right)}{\phi_{A}^{2}-\phi_{B}^{2}}.$
(8b)

We equate these expressions to eliminate $\chi$ and multiply through by a
factor $2(\phi_{A}^{2}-\phi_{B}^{2})$ to arrive at

$\begin{split}&\left[\frac{1}{N}\ln\left(\frac{\phi_{A}}{\phi_{B}}\right)-\ln\left(\frac{1-\phi_{A}}{1-\phi_{B}}\right)\right](\phi_{A}+\phi_{B})=\\\
&2\left(\frac{1}{N}-1\right)(\phi_{A}-\phi_{B})-2\ln\left(\frac{1-\phi_{A}}{1-\phi_{B}}\right).\end{split}$
(9)

Organizing the terms, we can write (9) as:

$\begin{split}&\frac{\left(\phi_{A}+\phi_{B}\right)}{N}\ln\left(\frac{\phi_{A}}{\phi_{B}}\right)+(2-\phi_{A}-\phi_{B})\ln\left(\frac{1-\phi_{A}}{1-\phi_{B}}\right)=\\\
&2\left(\frac{1}{N}-1\right)(\phi_{A}-\phi_{B}),\end{split}$ (10)

or

$\begin{split}&\frac{\left(\phi_{A}+\phi_{B}\right)}{N(\phi_{A}-\phi_{B})}\ln\left(\frac{\phi_{A}}{\phi_{B}}\right)+\frac{2-\phi_{A}-\phi_{B}}{\phi_{A}-\phi_{B}}\ln\left(\frac{1-\phi_{A}}{1-\phi_{B}}\right)=\\\
&2\left(\frac{1}{N}-1\right).\end{split}$ (11)

If we define variables $y$ and $z$ as:

$\begin{split}&y=\frac{\phi_{A}-\phi_{B}}{\phi_{A}+\phi_{B}}\\\
&z=\frac{\phi_{A}-\phi_{B}}{2-\phi_{A}-\phi_{B}}\end{split}$ (12)

equation (11) becomes:

$\frac{1}{Ny}\ln\left(\frac{1+y}{1-y}\right)-\frac{1}{z}\ln\left(\frac{1+z}{1-z}\right)=2\left(\frac{1}{N}-1\right).$
(13)

The variables $y$ and $z$ are relative measures of the partitioning of polymer
and solvent, respectively, between phases, and they are already commonly used
as an order parameter that defines the phase composition differences going
back at least to Cahn and Hilliard Cahn and Hilliard (1958). The variable $z$
is similar to $y$ in that it contains the difference over the sum of solvent
in each phase,

$z=\frac{\phi_{sB}-\phi_{sA}}{\phi_{sA}+\phi_{sB}},$ (14)

where $\phi_{s}=1-\phi$ is the solvent volume fraction by incompressibility.
It should be noted, however, that these mappings to $y$ and $z$ are not
conformal. Nevertheless, since we have one degree of freedom, we will specify
the variable $y$, and then solve the above equation for $z$.

In this pursuit, we first can recognize the $\tanh^{-1}$ function, defined by
the identity

$\tanh^{-1}(y)=\frac{1}{2}\ln\left(\frac{1+y}{1-y}\right),$ (15)

to rewrite equation (13) as

$\frac{\tanh^{-1}(y)}{Ny}-\frac{\tanh^{-1}(z)}{z}=\frac{1}{N}-1.$ (16)

Figure 2: Polymer-solvent solution. Phase diagram for a single polymer-solvent
solution, comparing different degrees of polymerization ($N=10,50,100,500$).
(a) Shows the analytical function $z(y)$ that describes the solvent
partitioning as a function of the polymer partitioning between phases, (b)
shows the analytical function $\chi(y)$, and (c) maps the coexistence curve to
$\chi(\phi)$. The numerical solution overlaps perfectly with the exact
analytical solution derived here. For large $N$, the asymptotic approximation
first used by Flory Flory (1942) matches most of the curve far from the
critical point. The red circles mark the critical points at each value of $N$.

If we define the function $h(x)$ as

$h(x)=\frac{\tanh^{-1}(x)}{x},$ (17)

then equation (16) becomes:

$h(z)=1+(h(y)-1)/N.$ (18)

We call this function, $h()$, the “FH function” since it features prominently
in our solution method for the Flory-Huggins model. Using the inverse FH
function, $h^{-1}(x)$, which returns the positive branch, we can solve the
equation for $z$ in terms of $y$:

$z=h^{-1}\left(\frac{1}{N}h(y)-\frac{1}{N}+1\right).$ (19)

Therefore, at this point, we have expressed the solution for the phase diagram
exactly in terms of $y$ and $z$. The function $h^{-1}(x)$ can be evaluated
easily using predefined lookup tables (our strategy), though other strategies
based on convergent series representations or asymptotic expressions of
$h^{-1}(x)$ are possible. We leave these approaches to future practitioners if
they are needed.

From the definition of $y$ and $z$, we can find $\phi_{A}+\phi_{B}$:

$\phi_{A}+\phi_{B}=\frac{2z}{y+z}.$ (20)

Finally, we can return to equation (8b), expressing $\chi$ in terms of $y$:

$\chi=\frac{\left(\frac{1}{N}-1\right)(\phi_{A}+\phi_{B})y-\ln\left(\frac{1-z}{1+z}\right)}{(\phi_{A}+\phi_{B})^{2}y},$
(21)

where we substitute in for $z$ and $\phi_{A}+\phi_{B}$ in terms of $y$ from
equations (19) and (20).

As mapped out in the schematic in Figure 1a, to construct the binodal curve,
we specify a value of $y$ between 0 and 1. Then, we find $z$ by applying
equation (19), then $\phi_{A}+\phi_{B}$ by applying equation (20), then the
corresponding $\chi$ value by applying equation (21). Note that we do not need
to explore negative values of $y$ since we arbitrarily assert that the $A$
phase is enriched in the polymer relative to the $B$ phase. If we wish to cast
the results in terms of $\phi_{A}$ and $\phi_{B}$ we can write:

$\begin{split}&\phi_{A}=\frac{z(1+y)}{z+y}\\\
&\phi_{B}=\frac{z(1-y)}{z+y}.\end{split}$ (22)

If we enumerate multiple values of $y$ between 0 and 1, we can construct the
full phase diagram. Our solution is therefore an implicit function, $\chi(y)$,
which we can map to $\chi(\phi)$.

The exact solutions to the Flory-Huggins theory are plotted in Figure 2 for
$N=$ 10, 50 , 100, and 500. In Figure 2(a-b), we show the functions $z(y)$ and
$\chi(y)$ for each value of $N$, which we then map to the function
$\chi(\phi)$ in Figure 2c. We confirm that the exact analytical solutions
match numerical solutions of the coupled nonlinear equations. Surprisingly,
the approximation first used by Flory Flory (1942) is able to capture the
asymptotic behavior of the phase distribution, especially for large $N$. It
seems that this approximation has not received significant attention, since
most attention has been paid to analytical approximations in the region near
the critical point and extensions from the critical point. For these reasons,
we briefly explain the approximation method of Flory, as it shares some
similarities to the implicit method we have applied and may be of practical
use.

###  II.3 Extending Flory’s approximation

Interestingly, in his 1942 article, Flory solved his theory in a rather
ingenious way Flory (1942) using an asymptotic arguments of the osmotic
pressure difference between solutions. Recognizing that the osmotic pressure
of the dilute phase would go to zero relative to the condensed phase, Flory
solved the osmotic pressure equation for $\chi$ keeping only terms containing
the condensed phase volume fraction. Specifically, he recognized that in the
limit $N\gg 1$, and away from the critical point, it is natural to expect that
$\phi_{B}\ll 1$, i.e., the solution is very dilute in polymer, and
$\phi_{A}\approx 1$, i.e., the second phase is rich in polymer. In this limit,
it follows that equation (7b) simplifies to:

$\left(\frac{1}{N}-1\right)\phi_{A}-\ln\left(1-\phi_{A}\right)-\chi\phi_{A}^{2}=0.$
(23)

In other words, there is an implicit relation giving $\phi_{A}(\chi,N)$ in the
form

$\chi=\frac{\left(\frac{1}{N}-1\right)\phi_{A}-\ln\left(1-\phi_{A}\right)}{\phi_{A}^{2}}.$
(24)

Flory used this relationship in order to approximate the $\chi$ that
corresponds to a given dense phase concentration ($\phi_{A}$), then
iteratively solved for the dilute phase concentration ($\phi_{B}$) using the
chemical potential constraint.

Interestingly, the starting approximation in equation (24) (with
$\phi_{A}\gg\phi_{B}$) can be applied to the chemical potential equation (7a)
directly to arrive at an implicit approximation for the dilute volume
fraction, $\phi_{B}$,

$\frac{1}{N}\ln\phi_{B}=\frac{1}{N}\ln\phi_{A}-\ln\left(1-\phi_{A}\right)-2\chi\phi_{A}$
(25)

or

$\phi_{B}=\frac{\phi_{A}}{\left(1-\phi_{A}\right)^{N}}e^{-2N\chi\phi_{A}},$
(26)

where $\chi$ for a specified $\phi_{A}$ is given in equation (24). This last
result illustrates that typically $\phi_{B}$ is exponentially small.

As seen in Figure 2c, this asymptotic formula works incredibly well to
describe the coexisting phases, especially for large $N$, with the exception
of a small region near the critical point. Note that the approximation is
computed for $2\phi_{c}<\phi_{A}<1$, with $\phi_{c}$ defined by equation 5.
Other approximations near the critical point are well documented, and they can
also be extended to regions further from the critical point Qian  _et al._
(2022); van Leuken  _et al._ (2023a).

Now that we have solved the single-polymer case, we can turn our attention to
what we can learn about multicomponent mixtures.

##  III Multicomponent mixtures

Most polymer solutions contain more than one polymeric species. Even for
solutions containing a polymer of one type, the polymer sample often contains
polymers of differing lengths (molecular weight), and each polymer of a
specified length must be treated as a separate species Huggins and Okamoto
(1967); van Leuken  _et al._ (2023b).

Here, we will introduce the extension of the Flory-Huggins model to multiple
components. Then, the procedure for implementing the implicit substitution
method for this system is expounded. Since the number of variables scales with
the number of components, inevitably, we must introduce numerous composite
variables along the way, which we define as needed. However, the solution
strategy remains the same—working with variables $z$ and $y_{i}$, then using
substitution of $\chi$ to arrive at one nonlinear equation with one unknown,
and then specifying only the allowed number of composite variables to arrive
at a valid solution.

###  III.1 Theoretical background

For a multicomponent system with $M$ polymeric species of volume fraction
$\phi_{i}$ and size $N_{i}$ in a solvent, the free energy density is Mao  _et
al._ (2019):

$\begin{split}\tilde{f}&=\sum_{i}\frac{\phi_{i}}{N_{i}}\ln(\phi_{i})+\left(1-\sum_{i}\phi_{i}\right)\ln\left(1-\sum_{i}\phi_{i}\right)\\\
&-\chi\sum_{i}\sum_{j}\alpha_{ij}\phi_{i}\phi_{j},\end{split}$ (27)

where the sums are over all $M$ polymeric components (not including the
solvent). The incompressibility constraint has been applied to eliminate the
solvent density, and linear terms in $\phi_{i}$ are not included, since they
will not influence thermodynamic equilibrium. The coefficients $\alpha_{ij}$
indicate the shape of the symmetric effective interaction matrix, while the
prefactor $\chi$ encodes the magnitude of interactions. We assume that the
matrix $\boldsymbol{\alpha}$ is known. Up to a linear term, the free energy
reduces down to the form in equation (1) if there is only one polymeric specie
with $\alpha_{11}=1$.

The chemical potential of each species,
$\mu_{i}=\partial\tilde{f}/\partial\phi_{i}$, is up to an additive constant:

$\mu_{i}=\frac{1}{N_{i}}\ln(\phi_{i})-\ln\left(1-\sum_{j}\phi_{j}\right)-\chi\sum_{j}2\alpha_{ij}\phi_{j}.$
(28)

The osmotic pressure of the solution is computed as
$\Pi=-\tilde{f}+\sum_{i}\phi_{i}\partial\tilde{f}/\partial\phi_{i}$:

$\Pi=\sum_{i}\left(\frac{1}{N_{i}}-1\right)\phi_{i}-\ln\left(1-\sum_{i}\phi_{i}\right)-\chi\sum_{i}\sum_{j}\alpha_{ij}\phi_{i}\phi_{j}.$
(29)

We will assume that two phases $A$ and $B$ in equilibrium have composition
$\mathbf{\Phi}_{A}$ and $\mathbf{\Phi}_{B}$, where $\mathbf{\Phi}$ indicates a
vector with the components consisting of the corresponding chemical volume
fractions. These phases satisfy the chemical equilibrium constraints of
$\mu_{iA}=\mu_{iB}$ and $\Pi_{A}=\Pi_{B}$. Here we will have $M+1$ equations
and $2M+1$ unknowns ($\mathbf{\Phi}_{A}$, $\mathbf{\Phi}_{B}$, and $\chi$).
Therefore, we can set only $M$ variables independently.

In order to make analytical progress simplifying the set of nonlinear
equilibrium constraints, we next cast the equations into a set of composite
variables. The following key composite variables are defined in the next
section: $y_{i}$, the partitioning of component $i$ in equation 30; $z$, the
partitioning of solvent in equation 31; $\beta_{i}$, the sum of the coexisting
volume fractions of component $i$ in equations 31 and 43 ; $\gamma_{i}$ the
differential of interactions of component $i$ between phases in equations 35
and 42; $\eta_{i}$ the differential of interactions of component $i$ between
phases relative to component 1 defined after equation 53, and finally $w_{i}$,
the relative partitioning of component $i$ relative to component 1 in equation
41. A reader who wishes to skip to the final working equations can skip to
equations 53 and 55, which together define the simplified master equation in
one variable, $z$.

###  III.2 Casting into composite variables

Following the pattern from the one-component case with the implicit
substitution method (Section II), we simplify the arguments of the $\ln$
functions by defining $y_{i}$ as

$y_{i}=\frac{\phi_{iA}-\phi_{iB}}{\phi_{iA}+\phi_{iB}},$ (30)

which are between -1 and 1; these variables represent the partitioning of
component $i$ between phases. Arbitrarily, we choose $0\leq y_{1}<1$, although
all other $y_{i}$ can take on negative or positive values.

The analogous variable for the solvent partitioning is denoted $z$:

$z=\frac{\sum_{i}\left(\phi_{iA}-\phi_{iB}\right)}{2-\sum_{i}\left(\phi_{iA}+\phi_{iB}\right)},$
(31)

which again corresponds to the difference divided by the sum of volume
fractions of solvent in respective phases, as in equation 14. If we define the
set of variables $\beta_{i}$,

$\beta_{i}=\phi_{iA}+\phi_{iB},$ (32)

as the sum of volume fraction of component $i$ in coexisting phases, then we
can relate $z$ and $y_{i}$:

$\sum_{i}y_{i}\beta_{i}=\left(2-\sum_{i}\beta_{i}\right)z.$ (33)

To start solving the system analytically, we first solve each constraint for
$\chi$. The $M$ chemical potential equations (28) give:

$\chi=\frac{\frac{1}{N_{i}}\ln\left(\frac{1+y_{i}}{1-y_{i}}\right)-\ln\left(\frac{1-z}{1+z}\right)}{\gamma_{i}}\quad\hbox{for
every $i$},$ (34)

with the expression in the denominator, $\gamma_{i}$, defined as:

$\gamma_{i}=\sum_{j}2\alpha_{ij}\beta_{j}y_{j}.$ (35)

Therefore, $\gamma_{i}$ is a measure of the differential of interactions of
component $i$ between phases.

From the osmotic pressure equation (29), we get:

$\chi=\frac{\sum_{i}\left(\frac{1}{N_{i}}-1\right)\beta_{i}y_{i}-\ln\left(\frac{1-z}{1+z}\right)}{\gamma_{0}}$
(36)

with $\gamma_{0}$ defined as:

$\begin{split}\gamma_{0}&=\sum_{i}\sum_{j}\alpha_{ij}\left(\phi_{iA}\phi_{jA}-\phi_{iB}\phi_{jB}\right)\\\
&=\sum_{i}\sum_{j}\alpha_{ij}y_{j}\beta_{i}\beta_{j}=\frac{1}{2}\sum_{i}\beta_{i}\gamma_{i}.\end{split}$
(37)

As in Section II, we can conveniently cast the equations for $\chi$ in terms
of $\tanh^{-1}$ functions:

$\begin{split}&\chi=\frac{\frac{2}{N_{i}}\tanh^{-1}(y_{i})+2\tanh^{-1}(z)}{\gamma_{i}}\quad\hbox{for
every $i$},\\\
&\chi=\frac{2\tanh^{-1}(z)+\sum_{i}\left(\frac{1}{N_{i}}-1\right)\beta_{i}y_{i}}{\gamma_{0}}.\end{split}$
(38)

For compact notation, we will define the variable $x$ as:

$x=\tanh^{-1}(z)$ (39)

so that $z=\tanh(x)$, and variables $x_{i}$ as

$x_{i}=\tanh^{-1}(y_{i}).$ (40)

Before proceeding we can simplify the expressions for $\gamma_{i}$ (equation
35) by defining a variable $w_{i}$, which is a measure of the relative
partitioning of component $i$ relative to the partitioning of component $1$,

$w_{i}=\frac{\beta_{i}y_{i}}{\beta_{1}y_{1}}=\frac{\phi_{iA}-\phi_{iB}}{\phi_{1A}-\phi_{1B}},$
(41)

with $w_{1}=1$. The definition of $\gamma_{i}$ in terms of the set $w_{i}$ is
therefore

$\begin{split}\gamma_{i}&=\beta_{1}y_{1}\sum_{j}2\alpha_{ij}w_{j},\quad\quad
i\neq 0\\\
\gamma_{0}&=\beta_{1}y_{1}\sum_{i}\sum_{j}\alpha_{ij}w_{j}\beta_{i}.\end{split}$
(42)

With the definition of $w_{i}$, by applying equation (33) we can solve for
$\beta_{i}$ in terms of the set of $y_{i}$, $z$, and $w_{i}$,

$\beta_{i}=\frac{2zw_{i}}{y_{i}\sum_{j}\left[w_{j}\left(1+\frac{z}{y_{j}}\right)\right]}.$
(43)

To summarize the mathematical manipulations thus far, we have recast the
equations in terms of composite variables that simplify the arguments of the
nonlinear logarithmic functions. Next, we work in terms of these composite
variables, until we can write one master equation with one unknown by
specifying a maximum of $M$ independent composite variables.

###  III.3 Deriving the master equation

The next step, using our implicit substitution method, is to equate the
equations for $\chi$. For this purpose, we first equate the $\chi$ from the
first component to all other components, to obtain $M-1$ identities. Then we
equate a convenient linear combination of $\chi$ from all components to the
$\chi$ specified by the osmotic pressure, in order to cancel the denominator
$\gamma_{0}$. Proceeding, the combined chemical potential equations give

$\frac{\frac{2}{N_{i}}x_{i}+2x}{\gamma_{i}}=\frac{\frac{2}{N_{1}}x_{1}+2x}{\gamma_{1}},$
(44)

which rearranged is

$x_{i}=x\left(\frac{\gamma_{i}}{\gamma_{1}}-1\right)N_{i}+\frac{\gamma_{i}}{\gamma_{1}}\frac{x_{1}N_{i}}{N_{1}}.$
(45)

We next apply the osmotic pressure constraint on $\chi$, which we write as:

$\frac{2x+\sum_{j}\left(\frac{1}{N_{j}}-1\right)\beta_{j}y_{j}}{\gamma_{0}}=\frac{\frac{2}{N_{i}}x_{i}+2x}{\gamma_{i}}.$
(46)

Multiplying this equation by $\beta_{i}\gamma_{i}$, then summing over all
components, we obtain

$\frac{\sum_{i}\beta_{i}\gamma_{i}}{\gamma_{0}}\left(2x+\sum_{j}\left(\frac{1}{N_{j}}-1\right)\beta_{j}y_{j}\right)=\sum_{i}\beta_{i}\left(\frac{2x_{i}}{N_{i}}+2x\right).$
(47)

Applying the definition of $\gamma_{0}$ so as to eliminate it, we find

$\left(2x+\sum_{j}\left(\frac{1}{N_{j}}-1\right)\beta_{j}y_{j}\right)=\sum_{i}\beta_{i}\left(\frac{x_{i}}{N_{i}}+x\right).$
(48)

Then, collecting terms in $x$, we get:

$\left(2-\sum_{i}\beta_{i}\right)x-\sum_{i}\beta_{i}y_{i}=\sum_{i}\frac{\beta_{i}x_{i}}{N_{i}}-\frac{\beta_{i}y_{i}}{N_{i}},$
(49)

which can be written compactly as:

$\left(2-\sum_{i}\beta_{i}\right)(x-z)=\sum_{i}\frac{\beta_{i}(x_{i}-y_{i})}{N_{i}}$
(50)

or as:

$\frac{\tanh^{-1}(z)}{z}=1+\frac{\sum_{i}\beta_{i}(\tanh^{-1}(y_{i})-y_{i})/N_{i}}{\sum_{i}\beta_{i}y_{i}},$
(51)

where we have reverted back to $y_{i}$ and $z$ from $x_{i}$ and $x$. Further,
substituting in for $\beta_{i}$ from equation (43) gives:

$\frac{\tanh^{-1}(z)}{z}=1+\frac{\sum_{i}w_{i}(\tanh^{-1}(y_{i})/y_{i}-1)/N_{i}}{\sum_{i}w_{i}}.$
(52)

We can again write this expression in terms of the FH function, $h()$, defined
by equation (17):

$h(z)=1+\frac{\sum_{i}w_{i}(h(y_{i})-1)/N_{i}}{\sum_{i}w_{i}}.$ (53)

Defining $\eta_{i}=\gamma_{i}/\gamma_{1}$, which is the differential measure
of interactions between phases for component $i$ relative to component 1, the
chemical potential constraints in equation (44) give:

$x_{i}=N_{i}(\eta_{i}-1)x+\eta_{i}N_{i}x_{1}/N_{1}$ (54)

or in terms of $y_{i}$:

$y_{i}=\tanh(N_{i}(\eta_{i}-1)\tanh^{-1}(z)+\eta_{i}N_{i}\tanh^{-1}(y_{1})/N_{1}).$
(55)

At this stage, we have reached our stated goal—to find one equation in one
unknown by specifying at most $M$ independent variables. With careful
inspection, one may see that specifying $y_{1}$ and the set $w_{i}$ for $i>1$
(which constitutes $M$ variables) is sufficient to achieve this objective.
Stated more mechanically, by substituting equation (55) into equation (53),
and specifying the $M$ composite variables $y_{1}$ and $w_{i}$ for $i>1$, we
have successfully transformed our system of $M+1$ equations into one master
equation, equation (53), in one unknown, $z$, the solvent partitioning between
phases.

However, this equation does not have a clear inversion formula, since $z$ is
contained differently as an argument of multiple nonlinear functions, not just
one nonlinear function as in the one polymer case. Therefore, in the most
general case, the master equation must be solved numerically, but there are
still significant advantages of this approach compared to the original
numerical challenge, as we characterize below. Note however, that there are
clear strategies to approximate the master equation solution near and far from
the critical point or when a particular pair interaction is dominant, as
outlined in Appendices A-D.

Figure 3: Polydisperse polymer-solvent solution. Comparison of the phase
diagram for two polydisperse samples composed of one polymer type. (a)
Coexistence curves plotting $\chi$ versus total volume fraction of polymer
$\phi=\sum_{i}\phi_{i}$. Each curve fixes the volume fraction of the condensed
phase, $\nu$. The binodal curves are calculated with a polymer with a uniform
molecular weight distribution between $N=10-100$ such that $M=91$, and
separately for a polymer with an exponential molecular weight distribution
between $N=10-100$, $\bar{\phi}_{i}/\bar{\phi}_{1}=\exp((N_{i}-10)/45)$. (b-c)
The molecular weight distributions of the overall solution (thick lines), and
the molecular weight distributions in the dense and dilute phases at
$\chi=0.8$. Note that for $\nu=0$ and $\nu=1$, the dilute and dense phases,
respectively, perfectly overlap with the overall molecular weight
distribution.

Next, instead of choosing $w_{i}$ directly, we choose the set of
$w_{i}y_{1}=(\phi_{iA}-\phi_{iB})/(\phi_{1A}+\phi_{1B})$ for $i>1$ and
$y_{1}$, so that the specified variables have the same denominator of
$\phi_{1A}+\phi_{1B}$, and so that they all go to zero at the global critical
point. From $w_{i}$, one can determine the set $\eta_{i}$. Then, the master
equation can be solved for $z$ numerically. With $z$ in hand, all values of
$y_{i}$ and $\beta_{i}$ can be calculated explicitly, then finally $\chi$ can
be calculated using one of the original equilibrium constraints.

Even though a numerical solution of the master equation is still needed in
general, a vast simplification has been achieved by all the algebraic
manipulations entailed by the implicit substitution method— turning the system
of $M+1$ nonlinear equations into one nonlinear master equation (53) with one
unknown. Again, we must work with the composite variables $y_{1}$ and $w_{i}$,
but as sketched in the schematic in Figure 1b-c, the solution can be directly
transformed back to $\mathbf{\Phi}$ space with the relation:

$\phi_{iA}=\frac{1}{2}\beta_{i}(1+y_{i})$ (56a)
$\phi_{iB}=\frac{1}{2}\beta_{i}(1-y_{i}).$ (56b)

Additionally, we must also check that the solutions returned are actually
physical, since admissible solutions to the master equation can return
unphysical negative values of the volume fractions.

###  III.4 Special cases

In what follows, we discuss applications of the multicomponent solution method
to three special cases. The first case is that of a single polymer with an
arbitrary number of different lengths, i.e., the case of polydispersity.
There, we find an exact analytical solution, and the master equation solution
can be completely cast in terms of two independent variables without any
numerical solution necessary. In the second case, we apply the solution method
to a two polymer, one solvent system, to demonstrate how the method can be
used to find the $\chi$ function over the full compositional space. Finally,
in the third case, we suggest strategies to tackle the problem of large
multicomponent systems of polymers using the implicit solution without
discretizing the entire composition space.

####  III.4.1 Polydisperse polymer of one type

If a solution consists of only a solvent and only one polymer of one type, but
that polymer is polydisperse with distinct values of $N_{i}$, then we can
often assume the interactions are equal for each length of polymer, i.e.,
$\alpha_{ij}=1$, which means that $\eta_{i}=\gamma_{i}/\gamma_{1}=1$ for all
$i$. In this case, the value of $y_{i}$ becomes independent of $z$:

$y_{i}=\tanh\left(\frac{N_{i}}{N_{1}}\tanh^{-1}(y_{1})\right).$ (57)

If this is the case, then we can again solve the equation for $z$ explicitly
with the inverse FH function $h^{-1}$, exactly as we did for the one-polymer
case:

$z=h^{-1}\left(1+\frac{\sum_{i}w_{i}(h(y_{i})-1)/N_{i}}{\sum_{i}w_{i}}\right).$
(58)

For arbitrary molecular weight distributions, the values of $w_{i}$ and
$y_{1}$ can be chosen independently to construct the full phase diagram.
However, in many cases, the molecular weight distribution of the polymer being
added to the solvent is known. If the average volume fraction of a particular
molecular weight in both phases, $\bar{\phi}_{i}$, can be written in terms of
$\nu$, the volume fraction of the phase $A$ relative to the total volume
Huggins and Okamoto (1967),

$\bar{\phi}_{i}=\nu\phi_{iA}+(1-\nu)\phi_{iB}=(\phi_{iA}-\phi_{iB})\left(\nu+\frac{1}{2y_{i}}-\frac{1}{2}\right),$
(59)

then the ratio between the average polymer density
$\bar{\phi}_{i}/\bar{\phi}_{1}$ is:

$\frac{\bar{\phi}_{i}}{\bar{\phi}_{1}}=\frac{\left(\nu+\frac{1}{2y_{i}}-\frac{1}{2}\right)}{\left(\nu+\frac{1}{2y_{1}}-\frac{1}{2}\right)}w_{i}.$
(60)

Solving for $w_{i}$, we get:

$w_{i}=\frac{\bar{\phi}_{i}}{\bar{\phi}_{1}}\frac{\left(\nu+\frac{1}{2y_{1}}-\frac{1}{2}\right)}{\left(\nu+\frac{1}{2y_{i}}-\frac{1}{2}\right)}.$
(61)

If $y_{1}$ and $\nu$ are specified, then the set $w_{i}$ is also fixed if the
overall molecular weight distribution is known. Therefore, we can reduce the
dimensionality of the phase diagram by only varying $y_{1}$ and $\nu$ in the
ranges $0<y_{1}<1$ and $0<\nu<1$, as schematically sketched out in Figure 1b.

In Figure 3, we present the phase diagram for two polymer samples of different
overall molecular weight distributions, one that is uniformly distributed
between $N=10-100$ and one that is exponentially distributed between
$N=10-100$, where for both samples the number of species is therefore $M=91$.
These curves are analytically calculated in Figure 3a, which is much better
than working with the 92 nonlinear chemical equilibrium equations we started
with. The curves are plotted as $\chi$ versus the overall polymer volume
fraction, $\phi=\sum_{i}\phi_{i}$, at fixed condensed phase volume fraction,
$\nu$. In Figure 3b-c, the molecular weight distributions of the samples and
the molecular weight distribution of the coexisting phases are plotted in the
form of a probability density function. When the dense polymer phase volume
fraction is small ($\nu\approx 0$), the binodal is pushed towards lower
overall polymer density, because the first appearing condensed phases
preferentially include polymers of higher molecular weight, so they dominate
the coexistence curve.

As the volume fraction of the condensed phase, $\nu$, increases, the left
branch of the coexistence curve (Figure 3) shifts to the right, since the
condensed phase is forced to more closely match the overall molecular weight
distribution as the condensed phase volume fraction increases. Although
$\nu\approx 1$ may be practically infeasible to reach, there, the minority
dilute phase must enrich polymers of low molecular weight since they are
heavily crowded in the dense phase. Note that at $\nu=0$ and $\nu=1$, the
molecular weight distribution of the dilute and dense phase, respectively,
perfectly match the molecular weight of the overall polymer sample. Finally,
comparing the two overall molecular weight distributions, the sample more
enriched in the longer polymers has a left branch of the coexistence curve
shifted to the left. In all cases, the critical point appears to be dominated
by the largest molecular weight polymer constituents, although it varies with
the specific details of the molecular weight distribution and value of $\nu$,
as previously found in Huggins and Okamoto (1967).

Also, similar approximations that reduce the compositional dimensionality can
be made when one pair interaction is assumed to be dominant, as explained in
detail in Appendices C and D, to arrive at another set of analytical
solutions.

Figure 4: A two-polymer-one-solvent system. Parts (a-c) correspond to the
phase diagram for a diagonally dominant interaction matrix, while parts (d-f)
show the same results for an off-diagonally dominant interaction matrix. (a)
The full phase diagram mapping the $\chi$ surface (with colormap corresponding
to the value of $\chi$) for all compositional space for fixed $\alpha_{ij}$ on
the $y_{1}$, $\tanh(y_{1}w_{2})$ plane. All regions connect at the critical
point at $y_{1}=0$ and $y_{1}w_{2}=0$, but are segmented by the singular
curves on the $\chi$ surface that are evident by thick black regions on the
$\chi$ surface. The three poles are computed on a small circle near the
critical point, and marked with white markers. These poles can be used to
demarcate the distinct differentiable regions of the $\chi$ surface. Note that
the colormap is saturated at the bounds of the color bar, but the contour
lines are linear between $\chi=-10$ and $\chi=10$ in increments of 0.1, though
negative values of $\chi$ do not appear in the feasible physical regions.
Black regions correspond to regions where no solution is available or where
$|\chi|>10$. (b) The regions of the $\chi$ surface which are physical, with
positive predicted polymer densities. (c) The binodal curve for fixed values
of $\chi$ transformed to $\mathbf{\Phi}$ space, either mapped out by the
gradient descent algorithm in Section IIIC or by extracting contours from the
physical regions of the $\chi$ surface in part 4b. The small gaps correspond
to the critical points of the binodal curves, which map to the same point in
the composite composition coordinates. (d-f) The same information as (a-c),
but for a different interaction matrix.

####  III.4.2 A two polymer, one solvent mixture

While we were able to find an exact analytical solution for an arbitrarily
polydisperse polymer, the two polymer-one solvent mixture has a master
equation that can only be solved numerically. Here the master equation
becomes:

$h(z)=1+\frac{(h(y_{1})-1)/N_{1}+w_{2}(h(y_{2})-1)/N_{2}}{1+w_{2}}$ (62)

with $y_{2}$ fixed as:

$y_{2}=\tanh(N_{2}(\eta_{2}-1)\tanh^{-1}(z)+\eta_{2}N_{2}\tanh^{-1}(y_{1})/N_{1})$
(63)

and $\eta_{2}$ defined as:

$\eta_{2}=\frac{\alpha_{12}+\alpha_{22}w_{2}}{\alpha_{11}+\alpha_{12}w_{2}}.$
(64)

To construct a phase diagram, we specify $y_{1}$ and $y_{1}w_{2}$, then solve
for $\eta_{2}$, then for $z$ solving the master equation numerically. From
there, we can specify $y_{2}$, $\chi$, and then solve explicitly for the
polymer densities in each phase. To discretize the phase diagram, we vary
$0<y_{1}<1$ and $-1<\tanh(y_{1}w_{2})<1$. Note that it is more regularly
behaved to plot the phase diagram in terms of $y_{1}$ and the product
$y_{1}w_{2}$, since they share the same denominator in terms of their
$\mathbf{\Phi}$ dependence.

In Figure 4, we present the solution for a two polymer-one solvent system for
two different shapes of the interaction matrix: (a-c) a diagonally dominant
interaction matrix and (d-f) an off diagonal dominant matrix. So as to not
push the coexisting densities to extreme values near 0 or 1, we choose
$N_{1}=4$ and $N_{2}=3$ as a starting point. The $\chi$ surface is calculated
over the full compositional landscape. Interestingly, all differentiable
regions emanate from the central critical point at $y_{1}=0$, $w_{2}y_{1}=0$,
but are separated by singular discontinuous curves in the phase diagram. Not
all of these regions or parts of these regions return physical, positive
densities, as shown in Figure 4b,e. Nevertheless, by transforming specified
contours with admissible solutions to $\mathbf{\Phi}$ space, we can map out
the phase coexistence curves in terms of $\phi_{1}$ and $\phi_{2}$, as
demonstrated in Figure 4c,f.

Clearly for $M=2$, if the value of the desired $\chi$ contour is known,
mapping out the full $\chi$ landscape may be more computationally expensive
than solving the equations in $\phi$ space for the fixed, desired value of
$\chi$. So the fully discretized $\chi$ surface implemented here are of
limited use for generating single binodal curves at fixed $\chi$.

A natural question at this point to ask is: how can one efficiently sample the
implicit $\chi$ surface to find a desired contour or sets of contours
corresponding to multiple coexisting phases? Further, how can one choose
sampling points so as to dodge the singular curves in the phase diagram?
Finally, can the sampling be efficiently scaled for large values of $M$? To
explore these questions, in the next section, we formulate such a sampling
strategy.

####  III.4.3 Finding contours in a high dimensional phase diagram

In the case of a high dimensional multicomponent solution, we have found an
implicit solution for $\chi$, requiring only one nonlinear solve in one
variable, $z$, to construct the surface:

$\chi=\chi(y_{1},\tanh(y_{1}w_{i}))=\chi(\mathbf{v}),$ (65)

where only the $M$ independent variables are listed as arguments in the vector
$\mathbf{v}$.

Without enumerating all possible inputs for the implicit function $\chi$, we
may not know beforehand where a given contour lies on the $\chi$ surface
corresponding to coexisting phases. If we wish to only find the contours of
this surface at fixed $\chi=\chi_{0}$ with minimal number of $\chi$
evaluations in a high dimensional phase diagram ($M\gg 1$), then we can employ
a gradient descent algorithm on the $\chi$ surface to arrive at the desired
contour. Defining the inputs of the $\chi$ function as a vector $\mathbf{v}$,
with step index $k$, the algorithm proceeds as:

$\mathbf{v}_{k+1}=\mathbf{v}_{k}-\lambda_{k}\nabla_{v}F(\mathbf{v}).$ (66)

The function $F$ is:

$F(\mathbf{v})=\left(\chi(\mathbf{v})-\chi_{0}\right)^{2},$ (67)

and its gradient can be computed by finite difference. The step size
$\lambda_{k}$ can be specified by a given algorithm, such as the Barzilai-
Borwein method Barzilai and Borwein (1988).

Once the desired contour is reached, one can then trace out a contour by
stepping along a direction, $\mathbf{u}$, perpendicular to the gradient
$\mathbf{g}=\nabla_{v}\chi$. Note, however, that there are $M-1$ such
coordinate directions that are perpendicular to $\mathbf{g}$. Here, we choose
just one of these coordinate directions,

$\mathbf{v}_{k+1}=\mathbf{v}_{k}-\lambda_{k}\mathbf{u}$ (68)

with

$\begin{split}&u_{i}=g_{i},\quad\quad i>1\\\
&u_{i}=-\frac{1}{g_{1}}\sum_{i>1}g_{i}^{2}\quad\quad i=1.\end{split}$ (69)

Here the step size $\lambda_{k}$ can be scaled by the magnitude of the
gradient vector or the location of the current iteration, $\mathbf{v}_{k}$.
The direction of the step away or towards the global critical point can be
controlled by the sign of $\lambda_{k}$. If ever the step size is too large
and the algorithm steps off of the fixed $\chi$ contour, the gradient descent
can be used to return to the contour.

In a high dimensional composition space, the binodal curves become higher
dimensional surfaces. Stepping in one tangential coordinate direction may not
sample the entire binodal surface on its own. Therefore, one would need to
choose a path that efficiently maps the boundary of the binodal surface, and
chooses the tangential directions based on the distance from the critical
point, the attributes of the gradient vector, and the historical location of
points already sampled. Here, we do not derive such a high dimensional search
algorithm.

The gradient descent algorithm works to converge to one contour in the
$\mathbf{v}$ space, corresponding to phase coexistence points in the
$\mathbf{\Phi}$ space. If the algorithm is initiated in a region that is
separated from the desired contour by a singular, non-differentiable curve in
the $\mathbf{v}$ space, then more than one starting point would be needed to
locate all matching contours in the phase diagram. If this is the case, then
it is natural to ask: what are efficient ways to choose candidate starting
points?

Notably, all regions connect at the global critical point at $\mathbf{v}=0$,
but the regions are separated by singularities which extend radially from the
critical point, slicing the phase diagram like a pie, as seen clearly in
Figures 4a and 4d. Therefore, the analytical approximations of the equations
by expanding near the critical point and finding poles can give a small number
of candidate starting points in each region located at the center of each pie
piece.

To sketch out the strategy to initiate a gradient descent algorithm, we next
turn to how to segment the phase diagram efficiently to generate candidate
starting points near the critical point.

Segmenting a hypersphere near the critical point:

Near the critical point, the partitioning of each species becomes negligible,
i.e., $|y_{i}|,|z|\ll 1$, and the master equation is:

$z^{2}=\frac{\sum_{i}w_{i}y_{i}^{2}/N_{i}}{\sum_{i}w_{i}}$ (70)

with the chemical potential constraint setting:

$y_{i}=N_{i}(\eta_{i}-1)z+\eta_{i}\frac{N_{i}}{N_{1}}y_{1}.$ (71)

Combined, these equations give a quadratic equation for $z$ which can be
solved exactly, as documented in Appendix B. Nevertheless, extra degrees of
freedom are still specified by the set $y_{1}$, $w_{i}$.

To eliminate these degrees of freedom, we propose the following strategy. We
look for candidate starting points on a small hypersphere $|\mathbf{v}|=v_{0}$
near the critical point, with $v_{0}\ll 1/\max(N_{i})$. First, we find the
poles of the $\chi$ function on this hypersphere. Then, we choose a small
number of points in regions demarcated by the poles of the $\chi$ function. By
choosing specified points between the poles, we guarantee that we are testing
each candidate region for the specified $\chi$ value.

For finite but small $z$ and $y_{i}$, the poles of the $\chi$ function are
defined by the equation:

$\gamma_{i}\propto\sum_{j}\alpha_{ij}w_{j}=0$ (72)

which corresponds to finding the null space of the $\mathbf{\alpha}$ matrix,
and then rescaling so that $w_{1}=1$. If $\mathbf{\alpha}$ is full rank, then
the null space is empty. However, one other denominator in the full
expressions can go to zero, causing singularities in the phase diagram,
namely:

$\sum_{j}\left[w_{j}\left(1+\frac{z}{y_{j}}\right)\right]=0.$ (73)

Near the critical point, these poles must be solved for numerically for the
values of $w_{i}$ and $y_{1}$ that lead to singularities, making use of the
exact solution for $z$ near the critical point.

Note that the poles can actually be interpreted as curves on a given
hypersphere. For $M=2$, they correspond to points on a circle. For $M=3$, they
correspond to arcs on a sphere, and so on. In order to segment the space in
higher dimensions, one would need to find the intersection of these singular
curves on the hypersphere near the critical point, to find the boundaries of
the nonsingular “continents” on the hypersphere.

Once we have found the poles near the critical point at fixed $v_{0}$, we can
segment the hypersphere radially between these poles, making sure to test at
least one point in the gradient descent algorithm in each “continent” of the
segmented hypersphere. This would ensure that our algorithm can find a desired
contour, if it exists, and could find multiple contours of the same function
that are separated by non-differentiable curves.

Applying the gradient descent algorithm:

Here, we explicitly apply the gradient descent algorithm only for $M=2$, to
find the contours in Figures 4c,f, validating the results by comparison to the
full $\chi$ surface. The poles of the function do segment the phase diagram
for $M=2$, as demonstrated by the white points marked in Figures 4a,d on a
circle a fixed distance from the critical point. Choosing a handful of
starting values between the singular points on the circle, we can trace out
the fixed values of $\chi$, and then transform to $\mathbf{\Phi}$ space, to
identify curves matching Figures 4c,f.

While we demonstrate the gradient descent algorithm for $M=2$, we leave the
exploration of $M>2$ to future work. Note that this gradient descent algorithm
is very advantageous for cases where $M>4$, where it may be difficult to
discretize the $\mathbf{v}$ space over the whole compositional range a priori
Mao  _et al._ (2019). In fact, while the number of function evaluations at
each point scales with $M$ (to evaluate the components of the gradient
vectors), and the number of candidate points may also increase with $M$, the
gradient descent algorithm is otherwise insensitive to $M$, which is much more
favorable than discretizing all combinatoric possibilities in composition
space. Still, the stepping along the binodal surfaces in high dimensions must
resolve enough points so as to map out the coexistance surface to significant
resolution, which may be a computationally demanding task.

##  IV Conclusions

Using the implicit substitution method, we are able to analytically define the
coexistence curves for a polymer-solvent system and a polydisperse polymer-
solvent system. For generalized mixtures of many components, we have
simplified a large system of nonlinear equations into one nonlinear equation
in one composite variable. We have presented methods to approximate the
solutions, not only for the single-polymer solvent system (Section IIC), but
also in the case of multicomponent solutions far from their critical point
(Appendix A) and close to their critical point (Appendix B). More analytical
solutions are possible in arbitrary mixtures if there is one dominant
interaction, where all other non-interacting species operate as non-
interacting crowders (Appendices C-D).

As a test of the multicomponent master equation, we have mapped out the
implicit $\chi$ surface for two polymeric components over the full
compositional space, and demonstrated the demarcation between differentiable
regions of the $\chi$ surface that defines the binodal curve. Further, we have
suggested strategies to find full binodal curves at fixed $\chi$ without
discretizing the full compositional space, which we believe will be very
useful for $M>4$, but we leave a full exploration of these methods to future
studies.

Ample opportunities exist to develop the methods presented here further.
Beyond improving the computational performance of the method for many
components, one could seek analytical perturbative series solutions to the
general master equation. Further, one could implement the method to large
multicomponent mixtures, and explore the shape of contours of the $\chi$
surface for different interaction matrices. Moreover, the approach can be
directly extended to more detailed equations of state beyond Flory-Huggins,
although compact analytical solutions of a final master equation are not
guaranteed for all equations of state. In Appendix E, we highlight possible
extensions of the implicit substitution method to other equations of state.

The composite composition variables are experimentally accessible if the
composition of the dilute and concentrated phases are simultaneously measured.
Therefore, experimental phase diagrams may also be mapped to composite
composition space for direct comparison to the theoretical $\chi$ surface for
a particular $\boldsymbol{\alpha}$ matrix at fixed temperature.

The implicit substitution method, while possibly un-intuitive at first sight,
has given us a new lens to explore the phase behavior of polymeric phase
separation. Here we demonstrated its power in simplifying systems of coupled
nonlinear thermodynamic equations that are often impossible to solve
explicitly.

##  V Acknowledgements

We acknowledge support from the Princeton Center for Complex Materials, an NSF
MRSEC (DMR-2011750), and NSF for grant DMS/NIGMS 2245850. J.P.D. is also
supported by the Omenn-Darling Princeton Bioengineering Institute – Innovators
(PBI2) Postdoctoral Fellowship. We would like to thank Hongbo Zhao, Tejas
Dethe, and Qiwei Yu for helpful discussions and comments.

##  Appendix A Approximations far from the critical point for multicomponent
mixtures

For multicomponent mixtures, it may be difficult to seek approximations, since
multiple assumptions of limiting values of densities may need to be made
simultaneously. While the limiting solution near the critical point is
analytical (by applying the quadratic formula in equation (70) as detailed in
Appendix B), one general approximation does not work in all regions far from
the critical point.

However, we can explore the limiting behavior when the major component, which
we arbitrarily specify as component 1, has a large excess in the condensed
phase compared to the dilute phase relative to other components. We also
choose a secondary component that exhibits the opposite partitioning compared
to component 1, but exhibits more partitioning compare to the other species.
Therefore, this boundary corresponds to the edge of the phase diagram in our
chosen coordinates, where $y_{1}\rightarrow 1$,
$w_{i}y_{1}\rightarrow\pm\infty$, or $z\rightarrow 1$.

In our first analysis, we start with the solvent as the secondary dominant
component. If this is the case, $y_{1}\rightarrow 1-\epsilon$, where
$\epsilon$ is vanishingly small, with a solvent dominated dilute phase,
meaning also $z\rightarrow 1$. We can effectively solve the system as if no
other polymers are present. In the limit of $\epsilon\rightarrow 0$, where
correspondingly $w_{i}$ for $i>1$ go to zero, we get

$h(z)\approx 1+{(h(y_{1})-1)/N_{1}},$ (74)

which has solutions in terms of the inverse FH function:

$z=h^{-1}\left(1+{(h(y_{1})-1)/N_{1}}\right).$ (75)

We can subsequently solve for the set $y_{i}$ describing the partitioning of
other polymers by applying equation (55). We can map this approximation to all
of composition space for a specific value of $y_{1}$, $w_{i}$, and then check
that the master equation is satisfied up to some tolerance to accept the
approximation prediction, or to reject it.

In many cases, though, the dilute phase in one polymer is the dense phase in
the other, and vice versa, for a given polymer pair. In other words, the
solvent may not be significantly partitioned between phases. If that is the
case, we can seek approximations by assuming that one of the polymers is the
“solvent” while the other is the polymer, effectively rescaling the theory in
these dominant composition coordinates. If we were to do this with component
$j$ as the solvent, and component $y_{i}$ the dominant partitioning component,
we would find:

$y_{j}=-h^{-1}\left(1+{(h(y_{i})-1)N_{j}/N_{i}}\right).$ (76)

where $y_{j}$ is chosen to have the opposite sign of $y_{i}$. All other
composite variables, $y_{k}$, can be computed by the analogue of equation (55)
that treats component $j$ as the “solvent” component that satisfies
incompressibility, namely:

$y_{k}=\tanh\left(-\frac{N_{k}}{N_{j}}\left(\frac{\eta_{k}}{\eta_{i}}-1\right)\tanh^{-1}(y_{j})+\frac{\eta_{k}N_{k}}{\eta_{i}N_{i}}\tanh^{-1}(y_{i})\right)$
(77)

where for the solvent $y_{0}=-z$.

We can use this procedure to map out the areas near the boundary of our
composition space, where $y_{1}\rightarrow 1$ or
$w_{i}y_{1}\rightarrow\pm\infty$. Further, by considering each possible pair
of polymers or solvent combination, we can generate a combinatoric set of
approximations that span all regions far from the critical point, that may at
least make for useful guesses for the full master equation. Similar to the
one-polymer case, we expect that these approximations will be most accurate
when there are significant differences in the polymer lengths.

##  Appendix B Approximations near the critical point

Near the critical point, the partitioning of each species becomes negligible,
i.e., $|y_{i}|,|z|\ll 1$, and the master equation is

$z^{2}=\frac{\sum_{i}w_{i}y_{i}^{2}/N_{i}}{\sum_{i}w_{i}},$ (78)

with the chemical potential constraint setting:

$y_{i}=N_{i}(\eta_{i}-1)z+\eta_{i}\frac{N_{i}}{N_{1}}y_{1}.$ (79)

This approximation is valid as long as:

$\begin{split}&z\ll 1\quad\&\quad z\ll\frac{1}{N_{i}(\eta_{i}-1)}\\\ &y_{1}\ll
1\quad\&\quad y_{1}\ll\frac{N_{1}}{N_{i}\eta_{i}}.\end{split}$ (80)

Combined, these equations give a quadratic equation for $z$ which can be
solved exactly. The quadratic equation explicitly is:

$\begin{split}0=&z^{2}\left(-1+\frac{\sum_{i}w_{i}N_{i}(\eta_{i}-1)^{2}}{\sum_{i}w_{i}}\right)\\\
&+z\left(\frac{\sum_{i}2w_{i}N_{i}\eta_{i}(\eta_{i}-1)y_{1}/N_{1}}{\sum_{i}w_{i}}\right)\\\
&+\frac{\sum_{i}w_{i}N_{i}\eta_{i}^{2}y_{1}^{2}/N_{1}^{2}}{\sum_{i}w_{i}}.\end{split}$
(81)

Picking the set $w_{i}$ and $y_{1}$, we solve for $z$, choosing the root that
gives positive $\beta_{i}$ for all $i$, if in fact either root gives positive
$\beta_{i}$.

##  Appendix C Adding non-interacting “crowders”

To add to our list of analytical solutions, we can consider the case of non-
interacting crowders added to a single strongly interacting polymer in
solvent. Although this ideal non-interacting assumption may be rare in
reality, it may be a reasonable approximation if one polymer self-interaction
is much larger in magnitude than all other interactions. We assume that the
dominant polymer is $i=1$, where $\alpha_{11}=1$. If non-interacting crowders
are added with species index: $2<i<M$, then since their corresponding rows and
columns $\alpha_{ij}=0$, we have $\eta_{i}=0$ for these species. Therefore,
the $\boldsymbol{\alpha}$ matrix only has one entry.

For the non-interacting polymers, we have:

$y_{i}=-\tanh\left(N_{i}\tanh^{-1}(z)\right).$ (82)

The master equation can be written as:

$\begin{split}h(z)-1=&\frac{\left(h(y_{1})-1\right)/N_{1}}{\sum_{i}w_{i}}\\\
&+\frac{\sum_{i=2}^{M}w_{i}\left(\frac{N_{i}\tanh^{-1}(z)}{\tanh\left(N_{i}\tanh^{-1}(z)\right)}-1\right)/N_{i}}{\sum_{i}w_{i}}\end{split}$
(83)

While $z$ is embedded in multiple nonlinear functions, $y_{1}$ appears as only
one argument of a nonlinear function. Therefore, instead of specifying $y_{1}$
and solving for $z$, we can specify $z$ and solve for $y_{1}$:

$\begin{split}h(y_{1})=&1+N_{1}\sum_{i}w_{i}\left(h(z)-1\right)\\\
&+\sum_{i=2}^{M}w_{i}\left(1-\frac{N_{i}\tanh^{-1}(z)}{\tanh\left(N_{i}\tanh^{-1}(z)\right)}\right)\frac{N_{1}}{N_{i}}.\end{split}$
(84)

In this case, $y_{1}$ can explicitly be expressed in terms of the inverse FH
function.

Similar to the polydisperse problem, we may know the overall ratios of the
crowders being added to the solution relative to the total solvent being
added. We can express this ratio in terms of $z$ and $\nu$ (the volume
fraction of the dense phase A):

$\begin{split}\bar{\phi}_{i}&=\nu\phi_{iA}+(1-\nu)\phi_{iB}=(\phi_{iA}-\phi_{iB})\left(\nu+\frac{1}{2y_{i}}-\frac{1}{2}\right)\\\
&=(\phi_{iA}-\phi_{iB})\left(\nu-\frac{1}{2\tanh(N_{i}\tanh^{-1}(z))}-\frac{1}{2}\right)\end{split}$
(85)

then the ratio between the average polymer density of component $i$ and the
solvent $\bar{\phi}_{i}/(1-\sum_{j}\bar{\phi}_{j})$ is:

$\frac{\bar{\phi}_{i}}{1-\sum_{j}\bar{\phi}_{j}}=\frac{\left(\nu-\frac{1}{2\tanh(N_{i}\tanh^{-1}(z))}-\frac{1}{2}\right)}{\left(-\nu+\frac{1}{2z}+\frac{1}{2}\right)}\frac{w_{i}}{\sum_{j}w_{j}}.$
(86)

If the left hand side of the equation (the ratio the crowder to solvent), is
known for $i>1$, and $\nu$ and $z$ are specified, then the above equations
give $M-1$ linear equations for the $M-1$ values of $w_{i}$ for $i>1$.

Therefore, if $\nu$ and $z$ are specified with a known total composition of
the crowders and solvent, then the set $y_{1}$ and $w_{i}$ are all specified.
To construct a phase diagram with crowders present, we can simply vary $\nu$
and $z$ between 0 and 1 to construct the full $\chi$ surface, then transform
to $\mathbf{\Phi}$ space.

Note that we may also use this “non-interacting” construction to include
finite “compressibility” of our solution, by including “holes” on the lattice
as one of the non-interacting species, along with the solvent and polymer
species. However, this lattice based approximation is typically a poor
approximation for the compressibility of hard-sphere liquids Frydel and Levin
(2012), and therefore can likely only approximate the behavior of real
liquids.

##  Appendix D  One dominant polymer-polymer interaction

Analogous to the case of one dominant self-interacting polymer is the case of
one dominant polymer-polymer pair interaction, meaning an interaction matrix
that is zero everywhere except for two off-diagonal entries. If we are to
reframe the theory such that one of the polymers is made into the “solvent”
species, and the other is made into the “dominant” interacting component, then
the approximation becomes equivalent to the previous section, where all other
species are only non-interacting crowders that are non-participating observers
to the main interacting components. In these reframed variables, the
interaction matrix would have only one non-zero entry, $\alpha_{11}$.

Note that for the case where the “solvent” species has length $N_{0}>1$, we
can rewrite the master equation as

$h(z)=1+\frac{\sum_{i}w_{i}(h(y_{i})-1)N_{0}/N_{i}}{\sum_{i}w_{i}}$ (87)

and the chemical potential constraints give $y_{i}$:

$y_{i}=\tanh(N_{i}(\eta_{i}-1)\tanh^{-1}(z)/N_{0}+\eta_{i}N_{i}\tanh^{-1}(y_{1})/N_{1}).$
(88)

Next, we could pursue the same approximations as in the previous sections
where $\eta_{i}=0$ for $i>1$, then frame the phase diagram in terms of $\nu$
and $z$. We leave an exploration of these approximations to future work.

It should be noted that there is a subtle difference between the “dominant-
interaction-approximation” in this section as compared to the “far from
critical” approximation in Appendix A. The dominant-interaction approximation
is valid only when the set $\eta_{i}\ll 1$ for $i\neq 1$, regardless of the
geometric distance from the critical point. However, for sufficiently strong
partitioning of a weakly interacting component $k$, $|w_{k}|\gg 1$, the
condition on $\eta_{k}$ will eventually break. So, in other words, the
“dominant-interaction approximation” breaks very far from criticality for
weakly interacting components. The “far from critical” approximation, on the
other hand, is even valid when the interaction matrix has many non-zero
entries, though its validity is only maintained far from the critical point.

To summarize all the approximations explored, building from Flory’s original
approach for a polymer-solvent system, we have derived an analytical
approximation far from the critical point. Next, we have derived an
approximation valid in a small region near the critical point. To bridge these
approximations, we have explored regions of the phase diagram where a
particular polymer-polymer or polymer-solvent interaction is dominant. While
more intricate approximations may be sought for triplets of species, we
suspect that for most regions in the phase diagram, the phase composition is
dominated by the two main interacting accumulated/depleted components between
the phases, and if not, we are either close to the global critical point or
very far from the global critical point—where other approximations are more
readily accessible.

##  Appendix E Extensions to other equations of state

Here, we can explore how this formulation could be extended to other equations
of state based on a lattice, focusing only on the polymer-solvent system. In
other theories, one could replace the $-\chi\phi^{2}$ term in the free energy
density by:

$-\chi\tilde{f}_{\mathrm{int}}(\phi)$ (89)

which is equivalent to assuming a $\phi$-dependent $\chi$ parameter.

In this case, from the chemical potential, we arrive at:

$\chi=\frac{\frac{2}{N}\tanh^{-1}(y)+2\tanh^{-1}(z)}{\tilde{f}_{\mathrm{int}}^{\prime}(\phi_{A})-\tilde{f}_{\mathrm{int}}^{\prime}(\phi_{B})}$
(90)

And from the pressure, we have:

$\chi=\frac{\left(\frac{1}{N}-1\right)(\phi_{A}+\phi_{B})y+2\tanh^{-1}(z)}{\tilde{f}_{\mathrm{int}}(\phi_{B})-\tilde{f}_{\mathrm{int}}(\phi_{A})+\phi_{A}\tilde{f}_{\mathrm{int}}^{\prime}(\phi_{A})-\phi_{B}\tilde{f}_{\mathrm{int}}^{\prime}(\phi_{B})}$
(91)

Substituting in:

$\begin{split}&\phi_{A}=\frac{z(1+y)}{z+y}\\\
&\phi_{B}=\frac{z(1-y)}{z+y}\end{split}$ (92)

and equating the values of $\chi$, we generally have a nonlinear equation for
one unknown ($z$) in terms of one specified variable ($y$)—reducing our system
of two equations down to one equation. While not guaranteed, one may be able
to decouple the $z$ and $y$ dependence, defining a function similar to the FH
function $h()$ for a given functional interaction form. While this procedure
is marginally better than solving the full system of two equations for the
polymer-solvent system, it may be especially helpful for the multicomponent
case. We leave these generalizations to future work.

## References

  * Huggins (1941) M. L. Huggins, The Journal of Chemical Physics 9, 440 (1941).
  * Flory (1942) P. J. Flory, The Journal of Chemical Physics 10, 51 (1942).
  * Cantow (1967) M. J. Cantow,  _Polymer fractionation_ (Academic Press Inc, 1967).
  * Barton (2018) A. F. Barton,  _Handbook of Poylmer-Liquid Interaction Parameters and Solubility Parameters_ (Routledge, 2018).
  * Milczewska and Voelkel (2006) K. Milczewska and A. Voelkel, Journal of Polymer Science Part B: Polymer Physics 44, 1853 (2006).
  * Peinemann  _et al._ (2007) K.-V. Peinemann, V. Abetz, and P. F. Simon, Nature Materials 6, 992 (2007).
  * Angelopoulou  _et al._ (2023) P. P. Angelopoulou, L. T. Kearney, J. K. Keum, L. Collins, R. Kumar, G. Sakellariou, R. C. Advincula, J. W. Mays, and K. Hong, Journal of Materials Chemistry A 11, 9846 (2023).
  * Callaway  _et al._ (2018) C. P. Callaway, K. Hendrickson, N. Bond, S. M. Lee, P. Sood, and S. S. Jang, ChemPhysChem 19, 1655 (2018).
  * Brangwynne  _et al._ (2015) C. P. Brangwynne, P. Tompa, and R. V. Pappu, Nature Physics 11, 899 (2015).
  * Berry  _et al._ (2018) J. Berry, C. P. Brangwynne, and M. Haataja, Reports on Progress in Physics 81, 046601 (2018).
  * Pappu  _et al._ (2023) R. V. Pappu, S. R. Cohen, F. Dar, M. Farag, and M. Kar, Chemical Reviews 123, 8945 (2023).
  * Overbeek and Voorn (1957) J. T. Overbeek and M. J. Voorn, Journal of Cellular and Comparative Physiology 49, 7 (1957).
  * Simha and Somcynsky (1969) R. Simha and T. Somcynsky, Macromolecules 2, 342 (1969).
  * Scheutjens and Fleer (1980) J. Scheutjens and G. Fleer, The Journal of Physical Chemistry 84, 178 (1980).
  * Tanaka (1989) F. Tanaka, Macromolecules 22, 1988 (1989).
  * Tanaka and Stockmayer (1994) F. Tanaka and W. H. Stockmayer, Macromolecules 27, 3943 (1994).
  * Semenov and Rubinstein (1998) A. N. Semenov and M. Rubinstein, Macromolecules 31, 1373 (1998).
  * Ruzette and Mayes (2001) A.-V. G. Ruzette and A. M. Mayes, Macromolecules 34, 1894 (2001).
  * Wang (2008) Z.-G. Wang, The Journal of Physical Chemistry B 112, 16205 (2008).
  * Wang (2017) Z.-G. Wang, Macromolecules 50, 9073 (2017).
  * Wessén  _et al._ (2021) J. Wessén, T. Pal, S. Das, Y.-H. Lin, and H. S. Chan, The Journal of Physical Chemistry B 125, 4337 (2021).
  * Sanchez (1984) I. C. Sanchez, Macromolecules 17, 967 (1984).
  * Tsenoglou and Papaspyrides (2001) C. J. Tsenoglou and C. D. Papaspyrides, Polymer 42, 8069 (2001).
  * Scheinhardt-Engels  _et al._ (2004) S. Scheinhardt-Engels, F. Leermakers, and G. Fleer, Physical Review E 69, 021808 (2004).
  * Qian  _et al._ (2022) D. Qian, T. C. Michaels, and T. P. Knowles, The Journal of Physical Chemistry Letters 13, 7853 (2022).
  * van Leuken  _et al._ (2023a) S. H. van Leuken, R. A. van Benthem, R. Tuinier, and M. Vis, Macromolecular Theory and Simulations 32, 2300001 (2023a).
  * Mao  _et al._ (2019) S. Mao, D. Kuldinow, M. P. Haataja, and A. Košmrlj, Soft Matter 15, 1297 (2019).
  * Jacobs and Frenkel (2017) W. M. Jacobs and D. Frenkel, Biophysical journal 112, 683 (2017).
  * Jacobs (2023) W. M. Jacobs, Journal of Chemical Theory and Computation 19, 3429 (2023).
  * Huggins and Okamoto (1967) M. L. Huggins and H. Okamoto, in  _Polymer Fractionation_ (Elsevier, 1967) pp. 1–42.
  * Wühr  _et al._ (2014) M. Wühr, R. M. Freeman, M. Presler, M. E. Horb, L. Peshkin, S. P. Gygi, and M. W. Kirschner, Current biology 24, 1467 (2014).
  * Doi (2013) M. Doi,  _Soft Matter Physics_ (Oxford University Press, USA, 2013).
  * Cahn and Hilliard (1958) J. W. Cahn and J. E. Hilliard, The Journal of Chemical Physics 28, 258 (1958).
  * van Leuken  _et al._ (2023b) S. H. van Leuken, R. A. van Benthem, R. Tuinier, and M. Vis, Journal of Physics: Materials 7, 015005 (2023b).
  * Barzilai and Borwein (1988) J. Barzilai and J. M. Borwein, IMA journal of numerical analysis 8, 141 (1988).
  * Frydel and Levin (2012) D. Frydel and Y. Levin, The Journal of Chemical Physics 137 (2012).