# What is the best simulation approach for measuring local density
fluctuations near solvo/hydrophobes?

Nigel B. Wilding  Robert Evans  Francesco Turci  H.H. Wills Physics
Laboratory, University of Bristol, Royal Fort, Bristol BS8 1TL, U.K.

###### Abstract

Measurements of local density fluctuations are crucial to characterizing the
interfacial properties of equilibrium fluids. A specific case that has been
well-explored involves the heightened compressibility of water near
hydrophobic entities. Commonly, a spatial profile of local fluctuation
strength is constructed from measurements of the mean and variance of solvent
particle number fluctuations in a set of contiguous sub-volumes of the system
adjacent to the solvo/hydrophobe. An alternative measure proposed by Evans and
Stewart (J. Phys.: Condens. Matter 27 194111 (2015)) defines a local
compressibility profile in terms of the chemical potential derivative of the
spatial number density profile. Using Grand Canonical Monte Carlo simulation,
we compare and contrast the efficacy of these two approaches for a Lennard-
Jones solvent at spherical and planar solvophobic interfaces, and SPC/E water
at a hydrophobic spherical solute. Our principal findings are that: (i) the
local compressibility profile $\chi({\bf r})$ of Evans and Stewart is
considerably more sensitive to variations in the strength of local density
fluctuations than the spatial fluctuation profile ${\cal F}({\bf r})$ and can
resolve much more detailed structure; (ii) while the local compressibility
profile is essentially independent of the choice of spatial discretization
used to construct the profile, the spatial fluctuation profile exhibits strong
systematic dependence on the size of the subvolumes on which the profile is
defined. We clarify the origin and nature of this finite-size effect.

##  I Introduction

Fluids adsorbed at substrates or extended solutes typically display
distinctive local density profiles together with enhanced density
fluctuations. Perhaps the best-known examples arise in the context of the
continuous wetting and drying surface phase transitions that occur for a
planar substrate of macroscopic interfacial area. A continuous wetting
transition corresponds to adsorption from the bulk vapor whereby a film of
liquid intrudes between the vapor and a (sufficiently attractive) substrate
whose thickness – and hence the Gibbs excess adsorption – diverges
continuously as a thermodynamic control parameter is varied. The approach to
the transition is accompanied by growing density-density correlations parallel
to the substrate, arising from capillary wave-like fluctuations, whose
correlation length $\xi_{\parallel}$ diverges, see e.g. the reviews by
Dietrich[1] and Bonn et. al[2]. Drying is the analogue of wetting, but now the
bulk is liquid, and a film of vapor intrudes between it and the (repulsive or
very weakly attractive) substrate. Recent theoretical and numerical studies by
Evans et. al.[3, 4] have shown that for realistic choices of fluid-fluid and
substrate-fluid potentials a continuous drying transition is likely to occur
at bulk coexistence by decreasing the strength of the attractive interaction
between the substrate and the fluid. At a continuous drying transition, the
thickness of the film of vapor and the negative of the Gibbs excess adsorption
diverge, accompanied by growing density-density correlations characterized by
a diverging correlation length $\xi_{\parallel}$. Other situations where a
diverging parallel correlation length occurs are: i) near a pre-wetting
critical point, ii) a fluid confined between two identical parallel planar
substrates on approaching a capillary critical point [5] and ii) a fluid in an
asymmetric planar slit where one substrate prefers to wet and the other
prefers to dry, e.g. Stewart and Evans[6] and references therein. For all
these cases there is ample understanding of the nature of the surface
criticality to enable proper analysis of simulation results.

The situation is different for adsorption at solutes of finite extent; the
optimal way of measuring solvent density fluctuations is not immediately
obvious. A suitable measure should treat the full range of solute size and
shape; solutes might range from simple near spherical molecules to much larger
objects such as colloidal particles and include more complex entities such as
proteins or biological molecules. Often, the solvent considered is liquid
water, and then one is concerned with hydrophobic solutes. Various attempts
have been made to connect enhanced density fluctuations in water with the
degree of hydrophobicity and also with drying phenomena; the recent review of
Rego and Patel[7] provides many references. Of course, one cannot define a
parallel correlation length for density-density correlations near a general
solute. Rather, one requires a robust measure that describes how the strength
of density fluctuations varies with the distance from the solute, the size of
the solute, and the proximity of the thermodynamic state point of the solvent
to bulk coexistence.

Experience with planar substrates indicates that there is a natural measure.
This is the local compressibility $\chi({\bf r})$ defined as the derivative of
the equilibrium density profile $\rho({\bf r})$ with respect to (w.r.t) the
chemical potential $\mu$, ie. $\chi({\bf r})=\partial\rho({\bf
r})/\partial\mu|_{T}$, which in bulk, where the density is constant, is
proportional to the usual thermodynamic quantity, the isothermal
compressibility. In planar substrate geometry, the maximum of $\chi(z)$
diverges (essentially) in the same way as $\xi_{\parallel}^{2}$ on approaching
a continuous surface transition. Employing classical Density Functional Theory
(DFT) calculations for Lennard-Jones (LJ) fluids, Evans and Stewart[8]
demonstrated that $\chi(z)$ provides a valuable measure for characterizing
hydrophobicity, or more generally solvophobicity, at planar substrates.
Specifically, they showed that for substrates where Young’s contact angle
$\theta$ is very large, but $\theta<\pi$, the maximum in $\chi(z)$ is at least
an order of magnitude larger than the bulk value $\chi_{b}$ and occurs at
distances $z$ within one or two atomic diameters of the substrate. They argued
that $\chi(z)$ is a much sharper indicator of the degree of solvophobicity of
a substrate, as measured by the contact angle, than is the extent of density
depletion, as extracted from the density profile $\rho(z)$. Subsequent papers
reinforced these ideas using Grand Canonical Monte Carlo (GCMC) simulations
and DFT for LJ fluids[9, 3] and simulations of SPC/E water[10]. Later studies
by Eckert et.al.[11, 12] introduced a (closely related) local thermal
susceptibility $\chi_{T}(z)=\partial\rho(z)/\partial T|_{\mu}$ that is the
temperature derivative of the density profile. Coe et. al.[13] showed that on
approaching critical drying, the singular behavior of $\chi(z)$ drives
identical singular behavior in $\chi_{T}(z)$.

Note that in an early paper [6], and references therein, $\chi(z)$ was termed
the local susceptibility since it is the direct analogue for continuum fluids
of the layer magnetic susceptibility $\chi_{n}$ in an Ising lattice subject to
a surface magnetic field. Specifically, $\chi_{n}=\partial m_{n}/\partial h$
is the derivative of the (average) magnetization $m_{n}$, in the $n$-th layer
away from the surface, w.r.t. the external magnetic field $h$; recall that $h$
is equivalent to the chemical potential $\mu$. The layer susceptibility
provides a powerful measure of the strength of local magnetization
fluctuations in layer $n$ and measurements of $\chi_{n}$, using Monte Carlo
simulations, have played a key role in elucidating the fundamental physics of
wetting transitions and phase transitions arising from confinement[14, 15].

Arguments that $\chi(r)$ provides an effective measure of the degree of
solvo/hydrophobicity at finite spherical solutes were presented in two recent
papers[16, 17] which draw upon ideas from earlier studies of critical drying
at (very large) spherical particles[18]. Similarly to planar systems, a
pronounced peak develops in $\chi(r)$ close to the solute whose height
increases with the proximity to bulk coexistence and now also with the radius
of the solute; the latter acts as a further variable in a comprehensive
scaling analysis of solvent thermodynamics and density fluctuations at a
solvophobe. Importantly, these studies, together with those for the planar
cases, emphasize the lack of water-specific mechanisms for the behavior near
an extended hydrophobe: the physics of the density fluctuations around a
generic solvophobe should be the same.

Of course, there are other measures of density fluctuations near hydrophobes
and Evans and Stewart[8] provide a summary. For models of water, or any fluid,
confined in a slit pore, the mean square fluctuation (variance) of the total
particle number provides a valuable measure of the overall compressibility of
the system and is often investigated in GCE simulation, see the commentary by
Bratko[19] and references therein. Here we focus on local measures. Acharya
et. al. [20] (see also Sarupria and Garde[21] and the review by Jamadagni et.
al. [22]) attempted a definition of a local compressibility that involves a
derivative of the density profile $\rho(z)$ w.r.t. pressure. However, it is
not clear precisely what the pressure is and how the derivative is performed.

In their paper, Acharya et. al. [20] also describe another quantity, denoted
as $\chi_{fl}(z)$, which measures the variance of particle number fluctuations
in a slab of certain thickness $\Delta z$ located at a distance $z$ from a
substrate. This quantity provides a simple means of defining a spatial
fluctuation profile and is essentially the quantity ${\cal F}(z)$ that we
shall define in Sec. II.1. Using molecular dynamics simulations, Willard and
Chandler[23] measured ${\cal F}(z)$ for SPC/E water near a planar Lennard-
Jones $12$-$6$ wall and found that this quantity increased, for $z$ close to
the wall, with decreasing wall-fluid attraction, i.e. with increasing contact
angle. They argued that ${\cal F}(z)$ provides a quantitative measure of the
degree of hydrophobicity. Evans and Stewart[8] noted that their DFT results
for $\chi(z)$ in the LJ fluid exhibited similar trends to those displayed in
the SPC/E results for ${\cal F}(z)$, but could not perform direct comparisons.

In the present work, we perform a systematic comparison of the utility of the
spatial fluctuation profile ${\cal F}({\bf r})$ (i.e. the generalization to
arbitrary geometry of ${\cal F}(z)$), and the local compressibility profile
$\chi({\bf r})$, for probing local density fluctuations. Our paper is arranged
as follows. Sec. II describes some pertinent background to fluctuations in
particle number. We consider the definitions of ${\cal F}({\bf r})$ and
$\chi({\bf r})$, emphasizing the distinction between these; $\chi({\bf r})$ is
proportional to the correlator (covariance) of the local particle number
$N({\bf r})$ and the total particle number $N$ whereas ${\cal F}({\bf r})$
measures the correlator of $N({\bf r})$ with itself. In Sec.III we present
GCMC results for three different physical systems: i) a LJ liquid confined
between two planar solvophobic walls, ii) a LJ liquid at a spherical
solvophobe, and iii) SPC/E water at a spherical hydrophobe. In each case, we
compare the relative merits of the two profile measures in quantifying the
strength and range of density fluctuations. We conclude in Sec.IV with a
summary and discussion.

##  II Background and methods

###  II.1 Spatial fluctuation profile

We start by considering the isothermal compressibility of a bulk (uniform)
fluid of volume $V=L^{d}$ (with $d$ the dimensionality of the system) defined
by $\kappa_{T}=-\frac{1}{V}\left(\frac{\partial V}{\partial p}\right)_{T}$,
with $p$ the system pressure. Within the grand canonical (constant $\mu VT$)
ensemble (GCE), $\kappa_{T}$ can be related to the configurational average
$\langle N\rangle$ and variance $\langle N^{2}\rangle-\langle N\rangle^{2}$ of
the fluctuating total number of particles $N$ by the well-known expression
[24]:

$\kappa_{T}=\frac{1}{k_{B}T\rho_{b}}\frac{\langle N^{2}\rangle-\langle
N\rangle^{2}}{\langle N\rangle}\>,$ (1)

where $\rho_{b}\equiv\langle\rho\rangle=\langle N\rangle/V$ is the fixed bulk
number density. Here we note that since $\langle N\rangle$ scales linearly
with $V$, then for $\kappa_{T}$ to be intensive requires that $\langle
N^{2}\rangle-\langle N\rangle^{2}$ scales similarly. This is equivalent to
requiring that the variance of the probability density function (pdf) of the
number density
$\langle\rho^{2}\rangle-\langle\rho\rangle^{2}=k_{B}T\rho_{b}^{2}\kappa_{T}/V$,
which is the standard result [25] for the finite-size scaling of Gaussian
density fluctuations that arise as a consequence of the central limit theorem.

For statistical ensembles in which the total particle number $N$ is fixed,
such as the constant-$NVT$ and constant-$NPT$ ensembles, it is clear that Eq.(
1) cannot be applied at the level of the total system. However, an estimate of
$\kappa_{T}$ can potentially be obtained by evaluating Eq. (1) over a domain
or ‘subvolume’[25]. In the limit in which the size of the subvolume $v$ is
large on the scale of the particle size, but small compared to the total size
of the system so that $v/V\to 0$, fluctuations in the number of subvolume
particles $N_{v}$ will effectively be grand canonical in form. This
observation has motivated several authors to attempt to evaluate via
constant-$N$ molecular dynamics simulations the spatial dependence of the
strength of density fluctuations in inhomogeneous fluids at substrates or
solutes by applying Eq. (1) to each of a set of contiguous subvolumes that
cover the spatial region of interest[23, 20, 21, 26]. Below we consider how
this strategy is implemented in practice for various solute/solvent geometries
of interest.

As a first example, consider the case of a three-dimensional ($d=3$) slit
geometry in which a solvent occupies the space between a pair of planar
substrates located at $z=0$ and $z=L$, with the system assumed periodic in the
$x$-$y$ plane so that the density profile varies only in the $z$ direction. It
is natural to discretize the system in the $z$ direction into identical thin
parallel slabs of thickness $\Delta z$, each having equal subvolume
$V(z)=L^{2}\Delta z$, where we have chosen the parallel area=$L^{2}$.
Evaluating the mean and variance of the solvent particle number fluctuations
in each slab yields a histogram of the fluctuation profile:

${\cal F}(z)=\frac{\langle N^{2}(z)\rangle-\langle N(z)\rangle^{2}}{\langle
N(z)\rangle}\>.$ (2)

Here we note in analogy to the discussion for $\kappa_{T}$ above, that since
$\langle N(z)\rangle\sim V(z)$, the numerator of Eq. 2 should similarly scale
linearly with $V(z)$ for ${\cal F}(z)$ to provide an intensive measure of the
local fluctuations.

Another common scenario considers an extended spherical solute particle fixed
at the origin and immersed in a solvent. For this geometry, it is natural to
discretize the space surrounding the solute into contiguous spherical shells
(concentric with the origin), each of which encompasses the space between some
radius $r$ and $r+\Delta r$ and has a subvolume $V(r)=4\pi r^{2}\Delta r$.
Evaluating the mean and variance of the solvent particle number fluctuations
within each shell yields a histogram of the radial fluctuation profile:

${\cal F}(r)=\frac{\langle N^{2}(r)\rangle-\langle N(r)\rangle^{2}}{\langle
N(r)\rangle}\>,$ (3)

for which similar scaling considerations as for $V(z)$ apply now with regard
to $V(r)$.

The approach can readily be extended to the case of an irregular solute, such
as a protein molecule, that lacks the symmetry of the above examples. Here,
the spatial variation of the fluctuation of the solvent density in $d=3$ can
be mapped using a fluctuation profile ${\cal F}({\bf r})$ calculated w.r.t. a
position vector ${\bf r}$ which we assume can range over the discrete lattice
vectors of a space-filling structure of eg. cubic subcells of volume $V({\bf
r})=(\Delta l)^{3}$. Evaluating $\langle N({\bf r})\rangle$ and $\langle
N^{2}({\bf r})\rangle$ for fluctuations in the solvent particle number in each
subcell yields a histogram ${\cal F}({\bf r})$:

${\cal F}({\bf r})=\frac{\langle N^{2}({\bf r})\rangle-\langle N({\bf
r})\rangle^{2}}{\langle N({\bf r})\rangle}\>.$ (4)

Each of the above approaches to defining a discretized fluctuation profile
entails a choice for the subvolume size and shape. The shape may be suggested
by the geometry of the problem, but the size needs to be sufficiently small to
resolve the pertinent features of the local density fluctuations. Furthermore,
when operating in a fixed $N$ ensemble, the subvolume should be much smaller
than the system size; otherwise, density fluctuations will be suppressed,
leading to biased results. However, it transpires that if one chooses a
subvolume that is smaller than some multiple of the diameter of the solvent
fluid particles, then the results can also be heavily biased. Aspects of the
issues of measuring fluctuations via subvolumes have previously been
highlighted by Villamaina and Trizac [27] and Román et al [28] who considered
the extent to which Eq. (1) when applied to a square subvolume of a two-
dimensional uniform fluid in the constant-$NVT$ ensemble yields an accurate
estimate for $\kappa_{T}$. It was found that it failed to do so for subvolumes
less than about $10$ particle diameters in linear size even in the limit when
this size was much smaller than that of the system as a whole [27]. Related
studies of uniform three-dimensional fluids have considered whether one can
correct measurements of subvolume compressibility to yield estimates of the
bulk compressibility [29, 30, 31].

The finding that subvolume estimates of $\kappa_{T}$ can exhibit serious
finite-size effects is a consequence of the central limit theorem: when the
subvolume size is insufficiently large, the local density fluctuations deviate
from the Gaussian form that pertains to very large subvolumes. This is true
even for situations that are far removed from a bulk or surface critical point
so that the correlation length for density-density fluctuations remains of
order the particle size. The consequence for fluctuation profiles is that the
numerator in each of Eqs. (2-4) scales nonlinearly with the subvolume, thus
engendering a subvolume size dependence of the fluctuation profile. However,
to date, the detailed consequences of this subvolume finite-size effect for
the sensitivity, resolution, and accuracy of measurements of local density
fluctuations in inhomogeneous fluids have not been investigated. We address
some of these here.

###  II.2 Local compressibility profile

The local compressibility profile as introduced by Evans and Stewart [8] is
defined within the GCE and takes the form

$\chi({\bf r})=\left.\frac{\partial\langle\rho({\bf
r})\rangle}{\partial\mu}\right|_{T}\>.$ (5)

Here $\langle\rho({\bf r})\rangle$ is the ensemble average of the
instantaneous density profile $\rho({\bf r})=\sum_{i=1}^{N}\delta({\bf r}-{\bf
r}_{i})$, where ${\bf r}$ is the d-dimensional position vector and $\mu$ is
the system chemical potential. It is straightforward to show (see appendix)
that for a bulk system for which $\chi({\bf r})=\chi_{b}$ is constant by
translational invariance, the local compressibility profile is related to the
bulk isothermal compressibility via $\chi_{b}=\rho_{b}^{2}\kappa_{T}$.

The definition eq. 5 is formally for a continuous $d$-dimensional density
profile. However, in practice $\chi({\bf r})$ is accumulated as a histogram by
discretizing the space of interest into a set of subvolumes of interest. Then
${\bf r}$ belongs to a discrete set of vectors running over the subvolumes,
which may, for example, be cubic cells, spherical shells, or planar slabs as
outlined in Sec. II.1. In the latter two cases, one obtains one-dimensional
profiles $\chi(r)$ and $\chi(z)$ respectively.

In simulations $\langle\rho({\bf r})\rangle$ is calculated as a
configurational average of histograms $\rho({\bf r})=N({\bf r})/V({\bf r})$
obtained by binning the solvent particle positions into the set of
subvolumes[10, 3] . A convenient approach for obtaining $\chi({\bf r})$ is to
explicitly perform the numerical derivative as per Eq. (5), by taking a finite
difference. The apparent need to perform two simulations at different $\mu$ to
achieve this can be neatly avoided by employing histogram reweighting [32] of
a single simulation. In short, a GCE simulation is performed at some chemical
potential $\mu$ of interest from which one accumulates an uncorrelated
sequence of $M$ simultaneous measurements of $\rho({\bf r})$ and the total
number of particles $N$. From this sequence, one first forms an estimate for
the ensemble-averaged density profile:

$\langle\rho({\bf r}|\mu)\rangle=\frac{1}{M}\sum_{i=1}^{M}\rho_{i}({\bf
r}|\mu)\;.$ (6)

Then this same sequence of measurements is reweighted w.r.t. a small notional
change $\Delta\mu$ in the chemical potential to yield an estimate for the
average profile corresponding to $\mu+\Delta\mu$:

$\langle\rho({\bf r}|\mu+\Delta\mu)\rangle=\frac{\sum_{i=1}^{M}\rho_{i}({\bf
r}|\mu)e^{\beta\Delta\mu N_{i}}}{\sum_{i=1}^{M}e^{\beta\Delta\mu N_{i}}}\>,$
(7)

where $N_{i}$ is the total number of particles for measurement $i$. Since
$\Delta\mu$ is notional, it can be chosen to be arbitrarily small (we use
$\Delta\mu=10^{-4}$). The local compressibility profile follows simply as

$\chi({\bf r}|\mu)=\frac{\langle\rho({\bf
r}|\mu+\Delta\mu)\rangle-\langle\rho({\bf r}|\mu)\rangle}{\Delta\mu}\>.$ (8)

A straightforward derivation starting from the grand partition function (see
the appendix) shows that $\chi({\bf r})$ can also be expressed in terms of a
correlator[10, 3] which takes the form

$\displaystyle k_{B}T\chi({\bf r})$ $\displaystyle=$
$\displaystyle\langle\rho({\bf r})N\rangle-\langle\rho({\bf r})\rangle\langle
N\rangle$ (9) $\displaystyle=$ $\displaystyle\frac{\left[\langle N({\bf
r})N\rangle-\langle N({\bf r})\rangle\langle N\rangle\right]}{V({\bf r})}\>.$
(10)

This expression provides an alternative route to that just described for
measuring $\chi({\bf r})$. Furthermore, it serves to expose the similarities
and differences between the two approaches for probing local fluctuations
represented by the quantities $\chi({\bf r})$ and ${\cal F}({\bf r})$.
Comparing Eqs. (10) and (4) shows that $\chi({\bf r})$ and ${\cal F}({\bf r})$
are fundamentally distinct in character. Specifically, $\chi({\bf r})$
correlates the instantaneous density profile $\rho({\bf r})$, or subvolume
particle number $N({\bf r})$, with the instantaneous total particle number
$N$, while ${\cal F}({\bf r})$ correlates the subvolume particle number with
itself. Thus, the numerical value of the product $\langle N({\bf r})N\rangle$
far exceeds that of $\langle N^{2}({\bf r})\rangle$. Consequently the central
limit theorem (and, by extension, the linear scaling with $V({\bf r})$ at
fixed $V$ of the numerator in Eq. 10) is satisfied to a much greater degree
for $\chi({\bf r})$ than is the case for ${\cal F}({\bf r})$, eq. 4. As we
shall see, this means that $\chi({\bf r})$ produces the bulk compressibility
value even for a minimal choice of the subvolume size on which it is defined,
in sharp contrast to ${\cal F}({\bf r})$.

We end this Section by noting that the integral of the difference ($\chi({\bf
r})-\chi_{b}$) over the volume available to the fluid measures the surface
excess compressibility $\chi_{e}$, which is the derivative of the Gibbs excess
adsorption w.r.t. chemical potential. $\chi_{e}$ is proportional to the
difference between the variance of the total number of particles in the
inhomogeneous fluid and the corresponding variance in the bulk fluid at the
same chemical potential. This quantity provides a powerful measure of the
integrated strength of density fluctuations as demonstrated explicitly for
planar systems[8].

##  III Monte Carlo simulation results

We have investigated the relative merits of the spatial fluctuation profile
${\cal F}$ and the local compressibility $\chi$ via Monte Carlo simulations of
three distinct physical setups:

  1. (i)

A Lennard-Jones solvent confined by a pair of solvophobic planar walls; the
profiles $\chi(z)$ and ${\cal F}(z)$ normal to the walls are studied using
planar slab subvolumes.

  2. (ii)

A Lennard-Jones solvent in contact with a solvophobic spherical solute; the
radial profiles $\chi(r)$ and ${\cal F}(r)$ are studied using spherical shell
subvolumes.

  3. (iii)

SPC/E water [33] in contact with a hydrophobic spherical solute; the profiles
$\chi({\bf r})$ and ${\cal F}({\bf r})$ are studied using cubic subvolumes.

In each instance, we work within the GCE [34], employing the venerable
algorithm of Metropolis, Rosenbluth, Rosenbluth, Teller, and Teller[35], which
is celebrated in this special issue.

###  III.1 Fluctuations in a slit with solvophobic planar walls investigated
with planar slab subvolumes

We construct a slit geometry from a cubic simulation box of volume $L^{3}$ by
placing planar walls at $z=0$ and $z=L$, with periodic boundary conditions in
the $x$ and $y$ directions. For the solvent, we employ a $12$-$6$ Lennard-
Jones (LJ) fluid that is truncated at $r_{c}=2.5\sigma$ (where $\sigma$ is the
LJ diameter) and left unshifted. The chemical potential $\mu$ is tuned to the
conditions of bulk liquid-vapor coexistence: the (reduced)
$\mu^{*}=-3.44610$[36] at a (reduced) temperature $T^{*}=1.0=0.842T_{c}$. The
fluid interacts with both walls through a wall-fluid potential that is
infinitely repulsive at $z=0$ and $z=L$ and has a long-range attraction of the
modified $9$-$3$ LJ form previously studied by Evans et al. [3]. By making the
dimensionless wall-fluid attraction very weak (we chose $\epsilon_{w}=0.01$,
in the notation of eq. 5 of ref.[3]) we render the wall strongly solvophobic,
giving rise to a state that is very close to critical drying [3]. For such a
state one expects substantial depletion of the density at each wall and
greatly enhanced density fluctuations. Note that in our slit setup, the liquid
is metastable w.r.t. capillary evaporation, but this does not prevent us from
studying the near-drying region at the walls (see the DFT study in Fig. 8 of
ref. [8] and Fig. 6 of ref. [3]).

We discretize the space between the walls into $L/\Delta z$ contiguous planar
slab subvolumes of prescribed equal thickness $\Delta z$ and volume
$V(z)=L^{2}\Delta z$. Histograms of the instantaneous density profile are
formed as $\rho(z)=N(z)/V(z)$. Fig. 1 shows our estimates of the ensemble
average $\langle\rho(z)\rangle$ (normalized by the independently determined
average bulk liquid density $\rho_{b}=0.654(1)$) for various choices of
$\Delta z$. The principal features of these profiles are a considerable
depletion in density near the solvophobic walls reflecting the incipient
drying regions, and a relaxation to the bulk liquid density far from the
walls. It should also be noted that within statistical uncertainties the
profiles are independent of $\Delta z$.

Figure 1: Monte Carlo results for the normalized density profile
$\langle\rho(z)\rangle/\rho_{b}$ calculated for a truncated LJ liquid at bulk
liquid-vapor coexistence and reduced temperature $T^{*}=1.0$ confined by a
pair of solvophobic planar walls. Results are shown for three choices of
subvolume size $V(z)=L^{2}\Delta z$ with $L=25\sigma$ and $\Delta z$ given in
the legend. The overall system volume is $V=25\sigma^{3}$ and the dashed red
line indicates $\langle\rho(z)\rangle=\rho_{b}$ where the bulk value is
determined from an independent simulation of a large fully periodic system.
Statistical errors are smaller than the symbol sizes.

We now present our measurements of the local compressibility profile $\chi(z)$
and the spatial fluctuation profile ${\cal F}(z)$. The definitions and
methodologies for calculating these quantities have been described in Sec. II.
Results for $\chi(z)$ are shown in Fig. 2(a) from which one sees (as already
established in ref. [3]) that $\chi(z)$ exhibits strong enhancement –in this
case by a factor of $~{}75$ compared to the bulk– in the region of the
incipient vapor-liquid interface that forms near the slit walls. Far from the
walls, $\chi(z)$ decays precisely to its bulk value
$\chi_{b}=\rho_{b}^{2}\kappa_{T}$ as denoted by the dashed red line
($\kappa_{T}=0.432(3)$ being determined by an independent simulation of the
liquid in a large periodic system). Similarly to $\langle\rho(z)\rangle$, we
see that the profiles for $\chi(z)$ are essentially independent of the choice
of $\Delta z$.

Matters are quite different for the spatial fluctuation profile ${\cal F}(z)$
shown in Fig. 2(b) and scaled by a factor of $(k_{B}T\rho_{b})^{-1}$ to allow
comparison with the bulk compressibility $\kappa_{T}$ whose value is indicated
by the red horizontal dashed line. ${\cal F}(z)$ was calculated using the same
sequence of configurations and slab subvolume sizes as for $\chi(z)$.
Nonetheless, the profile characteristics are quite distinct. Principally there
is a strong dependence on subvolume thickness $\Delta z$. Although all
profiles exhibit peaks near the wall, these are narrower than those of
$\chi(z)$ and grow with increasing $\Delta z$. The peak-to-trough enhancement
factor for ${\cal F}(z)$ is considerably less than for $\chi(z)$, showing that
the former is much less sensitive to variations in local density fluctuations
than the latter. It is also less accurate: the profile for
$(k_{B}T\rho_{b})^{-1}{\cal F}(z)$ fails to decay to its bulk value
$\kappa_{T}$ far from the walls. This latter observation accords with the bulk
studies of Villamaina and Trizac [27]. While further increasing the slab
subvolume thickness $\Delta z$ would help to ameliorate these issues with
${\cal F}(z)$, this would come at the cost of a reduction in the resolution of
the peak in this measure of density fluctuations.

Figure 2:  ${\bf(a)}$ Monte Carlo results for the local compressibility
profile $\chi(z)$ normalized by its bulk value $\chi_{b}$. These correspond to
the density profile shown in Figure 1. Results are shown for three choices of
subvolume size $V(z)=L^{2}\Delta z$ with $L=25\sigma$ and $\Delta z$ given in
the legend. The overall system volume is $V=25\sigma^{3}$ and the dashed red
horizontal line corresponds to $\chi(z)/\chi_{b}=1$. ${\bf(b)}$ The scaled
fluctuation profile $(k_{B}T\rho_{b})^{-1}{\cal F}(z)$ evaluated for the same
subvolumes considered in (a). Note the strong dependence on the subvolume
thickness $\Delta z$ and that the profile does not decay in the slit middle to
the correct value of the bulk compressibility $\kappa_{T}=0.432(3)$ in reduced
LJ units, as indicated by the dashed red horizontal line. In all cases,
statistical errors are smaller than the symbol sizes.

###  III.2 Fluctuations at a spherical solvophobic solute investigated with
spherical shell subvolumes

Next, we consider the case of a solvophobic spherical solute immersed in the
same LJ solvent used in sec. III.1. The center of the solute particle is fixed
at the origin and for the solute-solvent interaction we employ a potential of
the form given in Eq. (3) of Coe et al [17], with the attractive well depth
set to the very small (reduced) value $\epsilon_{sf}=0.01$ to render the
solute strongly solvophobic. We partition the volume around the solute into
concentric spherical shell subvolumes (centered on the solute) extending from
radius $r$ to $r+\Delta r$. The instantaneous radial density profile is
measured as $\rho(r)=N(r)/V(r)$ with $V(r)=4\pi r^{2}\Delta r$.

In Fig. (3) we show $\langle\rho(r)\rangle$ normalized by its bulk value
$\rho_{b}$ for two different solute particles of radii $R_{s}=3\sigma$ and
$R_{s}=5\sigma$. The profiles were accumulated using spherical shells of equal
thickness $\Delta r=0.0125\sigma$. One sees that in both cases close to the
solute particle there is a pronounced depletion in solvent density reflecting
the strong solvophobicity. The density profile in this region also exhibits
‘kinks’, particularly for the smaller solute. These occur on length scales of
about $\sigma$ and are remnants of the packing structure that occurs in bulk
fluids. Away from the solute, the density relaxes to its bulk value $\rho_{b}$
as indicated by the red dashed horizontal line.

Figure 3: Monte Carlo results for the radial density profiles for the
truncated LJ fluid at liquid-vapor coexistence for reduced temperature
$T^{*}=1.0$. The LJ solvent is in contact with a very weakly attractive
spherical solute and results are shown for solute radii $R_{s}=3.0\sigma$ and
$R_{s}=5.0\sigma$. The radial profiles were accumulated using spherical shell
subvolumes of thickness $\Delta r=0.125\sigma$. The overall system volume is
$V=25\sigma^{3}$ and the dashed red horizontal line indicates $\rho_{b}$, the
bulk value. Note the structure (‘kinks’) in the density profiles which stems
from packing effects. Statistical errors are smaller than the symbol sizes.

The corresponding local compressibility profiles $\chi(r)$ for the two solute
radii are shown in Fig. (4)a. Close to the solute there is a great enhancement
of the local compressibility compared to its bulk value. The peak in $\chi(r)$
for the larger solute particle extends over a greater range of radii than that
for the smaller one, in accord with previous findings [17]. Evident also is
considerable structure in the profile. Whilst this is inherited from the
density profile [3] the structure is much richer and more pronounced, with
subsidiary peaks in evidence. Far from the solute particle, the local
compressibility decays smoothly and accurately to its bulk value $\chi_{b}$.

The corresponding scaled spatial fluctuation profiles are shown in Fig. (4)b,
and were accumulated using the same subvolumes as for $\chi(r)$. Here one sees
that ${\cal F}(r)$ exhibits a much weaker response to the enhanced
fluctuations than does $\chi(r)$, and as in the slit case, the scaled profile
fails to decay to the value of the bulk compressibility. Here the failure is
dramatically clear.

Figure 4: (a) Monte Carlo results for the radial local compressibility profile
$\chi(r)$ for the truncated LJ fluid at liquid-vapor coexistence at
$T^{*}=1.0$ and in contact with a very weakly attractive spherical solute of
radius $R_{s}$. Results are shown for $R_{s}=3\sigma$ and $R_{s}=5\sigma$. The
subvolume shell thickness is $\Delta r=0.125\sigma$ and the dashed red
horizontal line corresponds to the bulk value. Note the structure in $\chi(r)$
which is related to the kinks in $\rho(r)$ seen in fig. 3. (b) Corresponding
results for the scaled radial fluctuation profile ${\cal F}(r)$ for the same
$\Delta r=0.125\sigma$. This profile resolves neither the extended range of
enhanced density fluctuations seen in $\chi(r)$ nor the detailed structure.
Moreover, it does not decay to the correct value of the bulk compressibility
$\kappa_{T}=0.432(3)$ in reduced LJ units, as indicated by the dashed red
horizontal line. In all cases, statistical errors are comparable with the
symbol sizes.

To investigate the effects of varying the subvolume shell thickness on the
results, we present in fig. 5 a comparison of profiles for $\chi(r)$ and
${\cal F}(r)$ for three different values of $\Delta r$. Fig. 5(a) shows that
$\chi(r)$ is insensitive to the subvolume size, and the detailed non-trivial
structure of $\chi(r)$ is reproduced consistently for each $\Delta r$. By
contrast the fluctuation profile ${\cal F}(r)$, Fig. 5(b) proves to be very
much less sensitive to the enhanced fluctuations near the solute surface than
$\chi(r)$. Specifically, ${\cal F}(r)$ has a much smaller peak compared to its
limiting value, fails to reflect accurately the range of the enhanced density
fluctuations, and does not resolve the structural features picked up by
$\chi(r)$. Whilst the sensitivity of ${\cal F}(r)$ increases with the
subvolume size, this is at the cost of spatial resolution and does not
approach the level of detail that is provided by $\chi(r)$.

Figure 5: (a) Monte Carlo results for the radial local compressibility profile
$\chi(r)$ for the truncated LJ fluid at liquid-vapor coexistence and
$T^{*}=1.0$ in contact with a very weakly attractive spherical solute of
radius $R_{s}=5.0\sigma$. The dashed red horizontal line corresponds to the
bulk value. Data are shown for three values of the spherical shell subvolume
thickness as indicated in the legend. ${\bf(b)}$ The scaled radial fluctuation
profile $(k_{B}T\rho_{b})^{-1}{\cal F}(r)$ evaluated for the same subvolumes
considered in (a). Note the strong dependence on the subvolume shell thickness
$\Delta r$ and that far from the solute the profile does not decay to the
correct value of the bulk compressibility $\kappa_{T}=0.432(3)$ (in reduced LJ
units), as indicated by the dashed red horizontal line. In all cases,
statistical errors are comparable with the symbol sizes.

Figure 6: Monte Carlo results for SPC/E water in contact with a very weakly
attractive spherical solute of radius $R_{s}=14$Å in a periodic box of volume
$V=(40$Å$)^{3}$. The temperature $T=300$K and $\mu$ is set to its bulk
coexistence value [37]. (a) The radial density profile
$\langle\rho(r)\rangle/\rho_{b}$; the inset shows a configurational snapshot.
(b) The local compressibility profile $\chi(r)$. (c) The scaled radial
fluctuation profile ${\cal F}(r)$ in units of ${\rm m}^{3}{\rm J}^{-1}$. Each
of $\langle\rho(\bf r)\rangle$, $\chi({\bf r})$ and ${\cal F}({\bf r})$ was
calculated for a 3D lattice of cubic subvolumes of size $\Delta
V=(0.5$Å$)^{3}$, and then averaged spherically to produce the radial profiles
shown. Statistical errors are comparable with the symbol sizes. The dashed red
horizontal line in (a) and (b) corresponds to the bulk value; that in (c)
corresponds to the bulk compressibility $\kappa_{T}=5.92\times 10^{-10}$
m3J-1. Note the lack of response of ${\cal F}(r)$ to enhanced density
fluctuations near the hydrophobic solute.

###  III.3 Fluctuations at a hydrophobic spherical solute investigated with
cubic subvolumes for SPC/E water

In this third example, we compare the local compressibility profile and the
spatial fluctuation profile in cubic subvolumes. This is a scenario that one
would likely adopt when seeking to map enhanced density fluctuations near an
irregularly shaped entity that exhibits hydrophobic regions, as is the case
for some large biomolecules.

The model that we have studied to illustrate this case is SPC/E water at a
spherical solute, which we have simulated in the GCE using the open-source
multipurpose Monte Carlo simulation engine DL_MONTE[38, 39]. The temperature
was set to $T=300$K and the chemical potential to its corresponding
coexistence value $\beta\mu=-15.24$ [37]. A solute of radius $R_{s}=14$Å was
fixed at the origin of a simulation box of size $V=(40$Å$)^{3}$ and the
solute-solvent interactions were assigned similarly to Sect. III.2 by
employing a potential of the form given in eq.(3) of Coe et al [17] between
the solute particle and the oxygen atom of the water molecules, with the
attractive well depth fixed, as in the previous example, by the very small
(reduced) value $\epsilon_{sf}=0.01$. The system was discretized into equal
cubic subcells of size $V({\bf r})=(0.5$Å$)^{3}$ and both the local
compressibility profile $\chi({\bf r})$ and the spatial fluctuation profile
${\cal F}({\bf r})$ were accumulated w.r.t. these. Although the subcells are
the fundamental units on which measurements were performed, we can improve
statistical sampling and generate a one-dimensional profile for comparison
purposes by exploiting the symmetry of our system to average both $\chi({\bf
r})$ and ${\cal F}({\bf r})$ spherically, yielding profiles for $\chi(r)$ and
${\cal F}(r)$. We stress, however, that the latter are not to be confused with
the profiles that would have resulted from choosing spherical shell subvolumes
in the manner of Sec. III.2.

The density profile for this system is shown in Fig. 6(a) and displays a $\sim
30\%$ depletion in the water density close to the solute and a weak
oscillatory decay to the bulk. The comparison of the spherically averaged
forms of $\chi({\bf r})$ and ${\cal F}({\bf r})$ is shown in Fig. 6(b) and (c)
respectively. In the former case, a strong enhancement of local
compressibility is seen, the spatial range of which correlates with the
depletion in the water density close to the solute, and which decays in a
similar oscillatory fashion to the bulk value far from the surface. By
contrast ${\cal F}({\bf r})$, while calculated in the same simulation and with
the same set of subcells as used to calculate $\chi({\bf r})$, shows no signal
of enhanced fluctuations at all –the spherically averaged profile
$(k_{B}T\rho_{b})^{-1}{\cal F}(r)$ rises from zero to a constant value which
is approximately $12$ times that of the bulk compressibility
$\kappa_{T}=5.92\times 10^{-10}$m3J-1. This complete insensitivity of ${\cal
F}({\bf r})$ to the density fluctuations is of course traceable to the small
cell size $V({\bf r})=(0.5$Å$)^{3}$ used in this case. While Fig. 6(b) shows
that this choice of cell size is necessary and warranted to resolve the
pertinent features of the local compressibility profile $\chi({\bf r})$, it
results in a complete lack of signal in ${\cal F}({\bf r})$.

One can better understand this finding by considering the statistics of very
small subvolumes[27]. When $V({\bf r})$ is smaller than the particle size,
cell occupancy is limited to $N({\bf r})=0$ or $1$. Accordingly, the
probability distribution function for occupation is binomial with some mean
occupancy $\langle N({\bf r})\rangle$. Given this, one readily finds that the
fluctuation profile eq. 4 evaluates as ${\cal F}({\bf r})=1-\langle N({\bf
r})\rangle$. When the subvolume cells are considerably smaller than the
particle size or the solvent density is low, as can be the case in proximity
to an extended solvo/hydrophobe, then $\langle N({\bf r})\rangle\ll 1$ and the
statistics approach the Poissonian limit for which the variance and mean of
the fluctuations in $N({\bf r})$ are equal. Accordingly, ${\cal F}({\bf r})$
approaches unity, which explains the absence of its response to density
fluctuations as reflected in Fig. 6(c). For position vectors ${\bf r}$
corresponding to bulk liquid SPC/E water, we find $\langle N({\bf
r})\rangle=0.0041(1)$ which is consistent with the limiting (large r) value of
$(k_{B}T\rho_{b})^{-1}{\cal F}({\bf r})=(k_{B}T\rho_{b})^{-1}(1-\langle N({\bf
r})\rangle)$ shown in Fig. 6(c) given the measured value of the bulk molecular
number density $\rho_{b}=3.348\times 10^{28}$m-3.

##  IV Summary and discussion

In this paper, we have compared the utility of two distinct approaches for
measuring the strength of local density fluctuations in fluids near
solvo/hydrophobes. The local compressibility $\chi({\bf r})$ introduced by
Evans and Stewart [8] in this context has been shown to provide high levels of
sensitivity, resolution, and accuracy regardless of the geometry or length
scales (subvolume size) for which it is applied. In particular as evidenced by
Figs. 4(a) and 6(b), $\chi({\bf r})$ can readily resolve the detailed features
of the local compressibility even when they occur on sub-particle length
scales. As mentioned in the Introduction, this measure is well rooted in the
statistical physics of interfacial phenomena, and its behavior in the vicinity
of surface phase transitions is well established.

By contrast, the spatial fluctuation profile ${\cal F}({\bf r})$, previously
employed by several authors to study density fluctuations in simulations of
water near hydrophobes, appears to be a considerably inferior measure in all
respects. Specifically, it generally exhibits a much weaker response to the
magnitude of local density fluctuations than $\chi({\bf r})$. This can only be
mitigated by increasing the subvolume size, which comes at the cost of a loss
of spatial resolution. The severity of this trade-off depends somewhat on the
geometry. For a planar solvo/hydrophobic substrate the subvolumes that one can
reasonably employ may not be so small as to render inadequate the sensitivity
and resolution of the spatial fluctuation profile. However, for spherical
solvo/hydrophobic solutes where the subvolumes in proximity to the surface are
necessarily small, the sensitivity and resolution are considerably worse than
those provided by the local compressibility. If one attempts to map the
spatial fluctuations around an irregular hydrophobe such as a bio-molecule
with acceptable resolution, then the spatial fluctuation profile will likely
fail to provide an adequate signal of the local density fluctuations, while
the local compressibility retains its full utility111Magnetic analogues of
$\chi(z)$ and ${\cal F}(z)$ have been defined by Binder and coworkers to
describe surface phase transitions in lattice-based magnetic models and termed
$\chi_{n}$ and $\chi_{nn}$ respectively. In accordance with our study, a much
stronger magnetic fluctuation signal occurs for $\chi_{n}$ than for
$\chi_{nn}$ [14]. However, in contrast to fluids, the choice of subvolume size
is uniquely predefined in the lattice context.

The differences between the two methods are traceable -in part at least- to
the degree to which they satisfy the central limit theorem in the limiting
case of a bulk fluid, as discussed in Sections II and III. On the smallest
length scales the occupation statistics of subvolumes are binomial, tending to
the Poissonian limit. Gaussian density fluctuations emerge in the bulk only
for sufficiently large subvolumes whose linear extent greatly exceeds the
local correlation length, and which contain correspondingly large numbers of
particles. Comparing the correlators that define the two approaches, eqs. (10)
and (4), one sees that the typical magnitude of the product $\langle N({\bf
r})N\rangle$ in the correlator for $\chi({\bf r})$ far exceeds that of
$\langle N^{2}({\bf r})\rangle$ appearing in ${\cal F}({\bf r})$. It follows
that the linear scaling, with subvolume size, of the numerator in the
correlator that is required to yield an intensive measure of fluctuation
strength is readily achieved for $\chi({\bf r})$ but not for ${\cal F}({\bf
r})$. This is why $\chi({\bf r})$ yields accurate estimates of the bulk
compressibility irrespective of the choice of subvolume size. And these
benefits of accuracy and high resolution extend beyond the case of purely bulk
fluctuations to yield greater sensitivity to the non-Gaussian near-critical
fluctuations that are the root cause of enhanced density fluctuations near
extended solvo/hydrophobes [17].

While the spatial fluctuation profile ${\cal F}({\bf r})$ has been employed in
several papers to study enhanced density fluctuation near extended hydrophobes
[23, 20, 21], we note that others [41, 22, 7] refer instead to the form of the
pdf of the subvolume particle number $P(N_{v})$. The appearance of enhanced
fluctuations in the form of a non-Gaussian tail at small $N_{v}$ has been
reported and attempts made to relate such behavior to the degree of
hydrophobicity of the substrate/ solute [41]. However, this tail is apparently
visible only for quite large subvolumes [7], which seems to rule out the use
of $P(N_{v})$ to create a high-resolution fluctuation profile.

Finally, we note that a key feature of the local compressibility $\chi({\bf
r})$ is that it is defined within the GCE. This open ensemble is doubtlessly
optimal for studying density fluctuations in inhomogeneous systems because
fluctuations can occur on all length scales up to and including the system
size. Furthermore, the GCE lends itself to accurate positioning of the
thermodynamic state of interest relative to (bulk) phase coexistence. Recall
that the deviation from bulk liquid-vapor coexistence, measured by the
chemical potential deviation, is a crucial ingredient in ascertaining the
origin, range and strength of solvo/hydrophobicity-induced density
fluctuations [3]. Additionally, as we have seen, the GCE offers practical
efficiencies such as the ability to calculate $\chi({\bf r})$ via histogram
reweighting (see Sec. II.2). Grand Canonical Monte Carlo is implemented in
several general-purpose molecular simulation engines such as DL_MONTE[39, 38]
and LAMMPS [42] and is thus readily accessible for a wide range of
applications including complex molecules. Nevertheless, it must be
acknowledged that many simulation studies of hydrophobic phenomena are not
performed grand canonically but rather using molecular dynamics in a closed
(constant-$N$) ensemble. While we expect the deficiencies of the spatial
fluctuation profile that we have identified in the GCE to be at least as
pronounced in closed ensembles as they are in the GCE, future work could
usefully consider the benefits offered by analogs of the local compressibility
to closed ensembles. Significant steps in this direction have recently been
reported by Eckert et al [12] who consider their local thermal susceptibility
in both the canonical and grand ensembles.

###### Acknowledgements.

NBW is grateful to Tom Underwood for guidance in the use of DL_MONTE. The
computer simulations were carried out using the computational facilities of
the Advanced Computing Research Centre, University of Bristol, as well as the
Isambard 2 UK National Tier-2 HPC Service (http://gw4.ac.uk/isambard/)
operated by GW4 and the UK Met Office, and funded by EPSRC (EP/T022078/1). RE
acknowledges support of the Leverhulme Trust, grant no.
EM-2020-029$\backslash$4

##  Appendix A Correlator for the local compressibility

For completeness, we include a derivation of Eq. 10 In the grand canonical
(constant-$\mu,V,T$) ensemble, the particle number $N$ fluctuates under the
control of the chemical potential $\mu$. Denote a configuration of $N$
particles as $\\{r_{i}\\}^{N}$ where $i=0\ldots N$. For an inhomogeneous fluid
under the influence of an external potential $V_{ext}({\bf r})$ that is a
function of the position vector ${\bf r}$, the Hamiltonian can be written
$H(\\{r_{i}\\}^{N}|\mu,V,T)=E(\\{r_{i}\\}^{N})+\int d{\bf r}\rho({\bf
r})V_{ext}({\bf r})-\mu N$, where $E(\\{r_{i}\\}^{N})$ is the configurational
energy. The grand partition function follows as
$\Xi(\mu,V,T)=\sum_{N=0}^{\infty}\int(\prod_{i=0}^{N}dr_{i})e^{-\beta H}$,
with $\beta=(k_{B}T)^{-1}$where we have suppressed phase space and
combinatorical factors.

The average density profile is obtained by a functional derivative of the
grand potential $\Omega=-\beta^{-1}\ln\Xi$ w.r.t. the external potential

$\displaystyle\frac{\delta\Omega}{\delta V_{ext}({\bf r})}$ $\displaystyle=$
$\displaystyle-\frac{1}{\beta\Xi}\frac{\delta\Xi}{\delta V_{ext}({\bf r})}$
(11) $\displaystyle=$
$\displaystyle\frac{1}{\beta\Xi}\sum_{N=0}^{\infty}\int(\prod_{i=0}^{N}dr_{i})\beta\rho({\bf
r})e^{-\beta H}$ (12) $\displaystyle=$ $\displaystyle\langle\rho({\bf
r})\rangle\>,$ (13)

which is a standard result.

Now define the local compressibility profile

$\displaystyle\chi({\bf r})$ $\displaystyle\equiv$
$\displaystyle\left.\frac{\partial\langle\rho({\bf
r})\rangle}{\partial\mu}\right|_{T}$ (14) $\displaystyle=$
$\displaystyle\sum_{N=0}^{\infty}\int(\prod_{i=0}^{N}dr_{i})\rho({\bf
r})\left[\frac{1}{\Xi}\beta Ne^{-\beta H}+e^{-\beta
H}\frac{1}{\Xi^{2}}\frac{\partial\Xi}{\partial\mu}\right].$ (15)

Thus

$\displaystyle k_{B}T\chi({\bf r})$ $\displaystyle=$
$\displaystyle\sum_{N=0}^{\infty}\int(\prod_{i=0}^{N}dr_{i})\left[N\rho({\bf
r})\frac{e^{-\beta H}}{\Xi}-\rho({\bf r})\frac{e^{-\beta H}}{\Xi}\langle
N\rangle\right]$ (17) $\displaystyle=$ $\displaystyle\langle\rho({\bf
r})N\rangle-\langle\rho({\bf r})\rangle\langle N\rangle\>.$ (18)

This result simplifies in the case of a bulk (uniform) system having constant
number density $\rho_{b}$, for which $\rho({\bf r})=\rho_{b}=\langle
N\rangle/V$ and $\chi({\bf r})=\chi_{b}$. Integrating eq. (18) w.r.t. the
differential volume $dV$, over the fixed system volume and noting that
$\int_{0}^{V}dV=V$ and $\int_{0}^{V}\rho({\bf r})dV=N$, one obtains

$\displaystyle Vk_{B}T\chi_{b}$ $\displaystyle=$ $\displaystyle\langle
N^{2}\rangle-\langle N\rangle^{2}$ (19)
$\displaystyle\frac{Vk_{B}T\chi_{b}}{\langle N\rangle}$ $\displaystyle=$
$\displaystyle\frac{\langle N^{2}\rangle-\langle N\rangle^{2}}{\langle
N\rangle}=\rho_{b}k_{B}T\kappa_{T}\>,$ (20)

where we have used Eq. (1). It follows that in the bulk

$\chi_{b}=\rho_{b}^{2}\kappa_{T}.$ (21)

## References

  * Dietrich [1988] S. Dietrich, “Wetting phenomena,” in  _Phase Transitions and Critical Phenomena_ , Vol. 12, edited by C. Domb and J. L. Lebowitz (Academic Press, 1988).
  * Bonn  _et al._ [2009] D. Bonn, J. Eggers, J. Indekeu, J. Meunier, and E. Rolley, “Wetting and spreading,” Rev. Mod. Phys. 81, 739–805 (2009).
  * Evans, Stewart, and Wilding [2017] R. Evans, M. C. Stewart, and N. B. Wilding, “Drying and wetting transitions of a Lennard-Jones fluid: Simulations and density functional theory,” J. Chem. Phys. 147, 044701 (2017).
  * Evans, Stewart, and Wilding [2019] R. Evans, M. C. Stewart, and N. B. Wilding, “From hydrophilic to superhydrophobic surfaces: a unified picture of the wetting and drying of liquids,” PNAS 116, 23901–8 (2019).
  * Evans and Parry [1990] R. Evans and A. O. Parry, “Liquids at interfaces: what can a theorist contribute?” J. Phys.: Condens. Matt. 2, SA15–32 (1990).
  * Stewart and Evans [2012] M. C. Stewart and R. Evans, “Phase behavior and structure of a fluid confined between competing (solvophobic and solvophilic) walls,” Phys. Rev. E 86, 031601 (2012).
  * Rego and Patel [2022] N. B. Rego and A. J. Patel, “Understanding hydrophobic effects: Insights from water density fluctuations,” Annual Review of Condensed Matter Physics 13, 303–324 (2022), https://doi.org/10.1146/annurev-conmatphys-040220-045516 .
  * Evans and Stewart [2015] R. Evans and M. C. Stewart, “The local compressibility of liquids near non-adsorbing substrates: a useful measure of solvophobicity and hydrophobicity?” J. Phys.:Condens. Matter 27, 194111 (2015).
  * Evans, Stewart, and Wilding [2016] R. Evans, M. C. Stewart, and N. B. Wilding, “Critical drying of liquids,” Phys. Rev. Lett. 117, 176102 (2016).
  * Evans and Wilding [2015] R. Evans and N. B. Wilding, “Quantifying density fluctuations in water at a hydrophobic surface: Evidence for critical drying,” PRL 115, 016103 (2015).
  * Eckert  _et al._ [2020] T. Eckert, N. C. X. Stuhlmüller, F. Sammüller, and M. Schmidt, “Fluctuation profiles in inhomogeneous fluids,” Phys. Rev. Lett. 125, 268004 (2020).
  * Eckert  _et al._ [2023] T. Eckert, N. C. X. Stuhlmüller, F. Sammüller, and M. Schmidt, “Local measures of fluctuations in inhomogeneous liquids: statistical mechanics and illustrative applications,” Journal of Physics: Condensed Matter 35, 425102 (2023).
  * Coe, Evans, and Wilding [2022a] M. K. Coe, R. Evans, and N. B. Wilding, “Measures of fluctuations for a liquid near critical drying,” Phys. Rev. E 105, 044801 (2022a).
  * Binder, Landau, and Ferrenberg [1995] K. Binder, D. P. Landau, and A. M. Ferrenberg, “Thin ising films with competing walls: A monte carlo study,” Phys. Rev. E 51, 2823–2838 (1995).
  * Binder, Landau, and Müller [2003] K. Binder, D. Landau, and M. Müller, “Monte carlo studies of wetting, interface localization and capillary condensation,” Journal of Statistical Physics 110, 1411–1514 (2003).
  * Coe, Evans, and Wilding [2022b] M. K. Coe, R. Evans, and N. B. Wilding, “Density depletion and enhanced fluctuations in water near hydrophobic solutes: Identifying the underlying physics,” Phys. Rev. Lett. 128, 045501 (2022b).
  * Coe, Evans, and Wilding [2023] M. K. Coe, R. Evans, and N. B. Wilding, “Understanding the physics of hydrophobic solvation,” J. Chem. Phys. 158, 034508 (2023), https://pubs.aip.org/aip/jcp/article-pdf/doi/10.1063/5.0134060/18140374/034508_1_5.0134060.pdf .
  * Stewart and Evans [2005] M. C. Stewart and R. Evans, “Critical drying at a spherical substrate,” J. Phys.: Condens. Matter 17, S3499–505 (2005).
  * Bratko [2007] D. Bratko, Faraday Discuss. 146, 382 (2010) 146, 382 (2007).
  * Acharya  _et al._ [2010] H. Acharya, S. Vembanur, S. N. Jamadagni, and S. Garde, “Mapping hydrophobicity at the nanoscale: Applications to heterogeneous surface and proteins,” Faraday Discuss. 146, 353–65 (2010).
  * Sarurpria and Garde [2009] S. Sarurpria and S. Garde, “Quantifying water density fluctuations and compressibility of hydration shells of hydrophobic solutes and proteins,” Phys. Rev. Lett. 103, 037803 (2009).
  * Jamadagni, Godawat, and Garde [2011] S. N. Jamadagni, R. Godawat, and S. Garde, “Hydrophobicity of proteins and interfaces: Insights from density fluctuations,” Annu. Rev. Chem. Biomol. Eng. 2, 147–71 (2011).
  * Willard and Chandler [2014] A. P. Willard and D. Chandler, “The molecular structure of the interface between water and a hydrophobic substrate is liquid-vapour like,” J. Chem. Phys. 141, 18C519 (2014).
  * Hansen and McDonald [2013] J. Hansen and I. R. McDonald,  _Theory of Simple Liquids_ , 4th ed. (Academic Press, 2013).
  * Rovere, Heermann, and Binder [1990] M. Rovere, D. W. Heermann, and K. Binder, “The gas-liquid transition of the two-dimensional Lennard-Jones fluid,” Journal of Physics: Condensed Matter 2, 7009 (1990).
  * Turci, Jack, and Wilding [2023] F. Turci, R. L. Jack, and N. B. Wilding, “Partial and complete wetting of droplets of active brownian particles,” (2023), arXiv:2310.07531 [cond-mat.stat-mech] .
  * Villamaina and Trizac [2014] D. Villamaina and E. Trizac, “Thinking outside the box: fluctuations and finite size effects,” European Journal of Physics 35, 035011 (2014).
  * Román, White, and Velasco [1998] F. L. Román, J. A. White, and S. Velasco, “Block analysis method in off-lattice fluids,” Europhysics Letters 42, 371 (1998).
  * Salacuse [2008] J. Salacuse, “Particle fluctuations within sub-regions of an n-particle, three-dimensional fluid: Finite-size effects and compressibility,” Physica A 387, 3073–3083 (2008).
  * Heidari  _et al._ [2018] M. Heidari, K. Kremer, R. Potestio, and R. Cortes-Huerto, “Fluctuations, finite-size effects and the thermodynamic limit in computer simulations: Revisiting the spatial block analysis method,” Entropy 20, 222 (2018).
  * Sevilla and Cortes-Huerto [2022] M. Sevilla and R. Cortes-Huerto, “Connecting density fluctuations and Kirkwood–Buff integrals for finite-size systems,” J. Chem. Phys. 156, 044502 (2022), https://aip.scitation.org/doi/pdf/10.1063/5.0076744.
  * Ferrenberg and Swendsen [1988] A. H. Ferrenberg and R. H. Swendsen, “New Monte Carlo technique for studying phase transitions,” Phys. Rev. Lett. 61, 2635–8 (1988).
  * Berendsen, Grigera, and Straatsma [1987] H. J. C. Berendsen, J. R. Grigera, and T. P. Straatsma, “The missing term in effective pair potentials,” J. Phys. Chem. 91, 6269–6271 (1987).
  * Frenkel and Smit [2002] D. Frenkel and B. Smit,  _Understanding Molecular Simulation_ , 2nd ed. (Academic Press, 2002).
  * Metropolis  _et al._ [2004] N. Metropolis, A. W. Rosenbluth, M. N. Rosenbluth, A. H. Teller, and E. Teller, “Equation of State Calculations by Fast Computing Machines,” J. Chem. Phys. 21, 1087–1092 (2004).
  * Wilding [1995] N. B. Wilding, “Critical-point and coexistence-curve properties of the Lennard-Jones fluid: A finite-size scaling study,” Phys. Rev. E 52, 602–611 (1995).
  * [37] “Liquid-vapor coexistence properties - SPC/E Water,” https://www.nist.gov/mml/csd/informatics/sat-tmmc-liquid-vapor-coexistence-properties-spce-water-lrc.
  * Brukhno  _et al._ [2021] A. V. Brukhno, J. Grant, T. L. Underwood, K. Stratford, S. C. Parker, J. A. Purton, and N. B. Wilding, “DL_MONTE: a multipurpose code for Monte Carlo simulation,” Mol. Sim. 47, 131–151 (2021).
  * [39] A. V. Brukhno and J. Grant and T. L. Underwood, “DL_MONTE source code v2.07,” {https://gitlab.com/dl_monte}.
  * Note [1] Magnetic analogues of $\chi(z)$ and ${\cal F}(z)$ have been defined by Binder and coworkers to describe surface phase transitions in lattice-based magnetic models and termed $\chi_{n}$ and $\chi_{nn}$ respectively. In accordance with our study, a much stronger magnetic fluctuation signal occurs for $\chi_{n}$ than for $\chi_{nn}$ [14]. However, in contrast to fluids, the choice of subvolume size is uniquely predefined in the lattice context.
  * Patel, Varilly, and Chandler [2010] A. J. Patel, P. Varilly, and D. Chandler, “Fluctuations of water near extended hydrophobic and hydrophilic surfaces,” J. Phys. Chem. B 114, 1632–7 (2010).
  * Thompson  _et al._ [2022] A. P. Thompson, H. M. Aktulga, R. Berger, D. S. Bolintineanu, W. M. Brown, P. S. Crozier, P. J. in ’t Veld, A. Kohlmeyer, S. G. Moore, T. D. Nguyen, R. Shan, M. J. Stevens, J. Tranchida, C. Trott, and S. J. Plimpton, “LAMMPS - a flexible simulation tool for particle-based materials modeling at the atomic, meso, and continuum scales,” Comp. Phys. Comm. 271, 108171 (2022).


# Modified S-wave $\pi\pi$ scattering amplitude for multiparticle PWA

I.A. Kachaev, VES experiment

(August 28, 2024)

###### Abstract

Suggested by Au, Morgan, Pennington (AMP) $S$-wave isospin $I=0$ $\pi\pi$,
$KK$ scattering amplitude is good enough to describe experimental data for the
moment. Still it has two disadvantages for use in multiparticle partial wave
analysis (PWA), namely sharp drop at the $KK$ threshold and unreasonable
behavior at $M(\pi\pi)>1.6\,GeV/c^{2}$. The drop is not seen in multiparticle
systems.

We suggest the modified AMP amplitude, mAMP, for the only aim, namely to
describe the broad part of $S$-wave $\pi\pi\,\to\,\pi\pi$ scattering in the
wide $M(\pi\pi)$ range in multiparticle PWA. The mAMP amplitude describes
threshold behavior of the $\pi\pi\,\to\,\pi\pi$ scattering and the wide
structure at $M\sim 1400\,MeV/c^{2}$ reasonably well. It is assumed that
narrow objects $f_{0}(980)$, $f_{0}(1500)$ are included in PWA separately. The
amplitude does not describe $\pi\pi\,\to\,KK$ scattering. The mAMP amplitude
is purely phenomenological.

(to be submitted to Physics of Atomic Nuclei)

##  1 Introduction

To perform multiparticle PWA we need to know two-particle amplitudes, the
_isobars_ in the PWA terminology. For $\pi\pi$ system in $P$, $D$, $F$ waves
two-particle amplitudes are mostly described by known resonances $\rho(770)$,
$f_{2}(1270)$, $\rho_{3}(1690)$ et al. Still in the $\pi\pi$ $S$-wave the
situation is more complicated. In addition to narrow objects $f_{0}(980)$,
$f_{0}(1500)$ there exists broad structure including at least threshold
peculiarities, wide $f_{0}(500)$, $f_{0}(1370)$ and probably some other
objects.

Scattering amplitude in the $\pi\pi$ $S$-wave was investigated in details in
AMP article [1]. This article is concentrated on the construction of the
scattering amplitude in the channels $\pi\pi$, $KK$ with proper analyticity
and unitarity. The best experimental data known at this time are used on
input. The amplitude is presented in the form appropriate for later use.

When AMP parametrization of $\pi\pi$ scattering was used in three-particle PWA
of $\pi^{-}\pi^{-}\pi^{+}$ system the following problems are found. The AMP
amplitude describes $\pi\pi$ $S$-wave in a whole, in the full range of
$M(\pi\pi)$. This leads to both physical and technical restrictions.

From the physical point of view amplitudes and phases of resonances in AMP
amplitude are fixed once and forever. The authors of [1] used so called
$P$-vector approach to describe production of the $\pi\pi$ system.

In the AMP amplitude there is a zero at $M(\pi\pi)\approx 1\,GeV/c^{2}$ due to
$f_{0}(980)$ object and the threshold in the $KK$ channel. If this zero is not
seen in the multiparticle system it should be compensated by the pole in
$P$-vector. It is impossible in the multiparticle PWA if in it a connection of
$\pi\pi$ system with production channel is described by a set of coupling
coefficients, as it is done in the Illinois PWA program [2].

From the technical point of view the AMP amplitude is parametrized via a set
of poles and background polynomial of 4th degree. It is known that polynomial
parametrization rapidly makes senseless outside of the range of definition.
Due to this fact AMP amplitude has unphysical maximum at $M(\pi\pi)\sim
1.6\,GeV/c^{2}$.

##  2 The $K$-matrix method

The AMP amplitude is constructed in $K$-matrix formalism [3, 4]. Formally the
transition from initial to final state $S_{ij}=\left<j|S|i\right>$ is
described by the unitary scattering operator, $SS^{+}=S^{+}S=I$.

Still it is much simple to work with Hermitian matrices then with unitary. Let
us construct the Hermitian matrix $K$ which describes $S$ exactly. After
separation of no-interaction part one can define the transition operator $T$
as

$S=I+2iT$ (1)

where $I$ is identity operator and factor $2i$ is introduced for convenience.
From unitarity of $S$ and the definition of $T$ we have

$T-T^{+}=2i\,T^{+}T=2i\,TT^{+}$ (2)

Introducing the inverse operator $T^{-1}$ we have

$(T^{+})^{-1}-T^{-1}=2i\,I\quad\mbox{so}\quad(T^{-1}+i\,I)^{+}=T^{-1}+i\,I$
(3)

Now we can introduce operators $K$ and $M$, from (3) they are Hermitian

$K^{-1}=T^{-1}+i\,I\quad\mbox{and}\quad M=K^{-1}$ (4)

It is known that if time-reversal invariance is hold than $K$ and $M$ matrices
are not only Hermitian but are also real and so symmetric. Explicit definition
for $T$ is

$T=K\left(I-iK\right)^{-1}=\left(I-iK\right)^{-1}K$ (5)

So defined transition amplitude $T$ is not Lorentz invariant. Let us construct
it in Lorentz invariant form. For decay of particle with mass $m$ to two
particles with masses $m_{a}$, $m_{b}$ phase space is $\rho=2q/m$ where $q$ is
breakup momentum

$\rho=\sqrt{\left[1-\left({\frac{m_{a}+m_{b}}{m}}\right)^{2}\right]\left[1-\left({\frac{m_{a}-m_{b}}{m}}\right)^{2}\right]}$
(6)

Phase space is normalized as $\rho\to 1$ at $m\to\infty$. By definition phase
space is considered a diagonal matrix, so in two channel case it is

$\rho=\left(\begin{array}[]{cc}\rho_{1}&0\\\ 0&\rho_{2}\end{array}\right)$ (7)

Lorentz invariant amplitude $\hat{T}$ is defined as

$T_{ij}=\\{\rho_{i}^{*}\\}^{1/2}\,{\hat{T}}_{ij}\,\\{\rho_{j}\\}^{1/2}$ (8)

or in matrix form

$T=\\{\rho^{+}\\}^{1/2}\,\hat{T}\,\\{\rho\\}^{1/2}$ (9)

Below the threshold of the channel its phase space become complex. Complex
conjugation in (8) is required for proper analytic continuation into the
region below the threshold of some channels.

Let us note a significant difference between $T$ and $\hat{T}$. In the one
channel case

$S=e^{2i\delta},\quad T=e^{i\delta}\sin\delta,\quad\hat{T}=\frac{1}{\rho}\,T$
(10)

where $\delta$ is phase shift. For any $\delta$ the amplitude $T$ is confined
in the circle with the centre $(0,i/2)$ and the radius $1/2$. The dependence
$T(s)$ is named the Argand diagram. The amplitude $\hat{T}$ has other
normalization and is used as Dalitz plot amplitude in PWA programs. On the
threshold of the reaction in the case of $S$-wave scattering $T\to 0$ but
$\hat{T}\to const$.

Now we can define Lorentz invariant matrix $\hat{K}$ as

$K=\\{\rho^{+}\\}^{1/2}\,\hat{K}\,\\{\rho\\}^{1/2}$ (11)

If $\hat{K}$ is taken to be real and symmetric, then $K$ is Hermitian and $S$
unitary even some $\rho$ becomes imaginary below the threshold of some
channels. From (4) we have

${\hat{K}}^{-1}={\hat{T}}^{-1}+i\rho$ (12)

Explicit form of $\hat{T}$ is

$\hat{T}=\hat{K}\left(I-i\rho\hat{K}\right)^{-1}=\left(I-i\hat{K}\rho\right)^{-1}\hat{K}$
(13)

##  3 The original AMP amplitude

In [1] (3.18) $\hat{K}$ matrix (named below $K$ for consistency with [1]) has
been parametrized as sum of poles and polynomial background (this is so called
$K$ solution)

$K_{ij}=\frac{(s-s_{0})}{4m^{2}_{K}}\sum_{p}\frac{f_{i}^{p}f_{j}^{p}}{(s_{p}-s)(s_{p}-s_{0})}+\sum_{n=0}c_{ij}^{n}\left(\frac{s}{4m^{2}_{K}}-1\right)^{n}$
(14)

By definition $m_{K}=1/2(m_{K+}+m_{K0})$. Here $s_{0}$ is known Adler zero
near the threshold of $\pi\pi$ system. We have found that amplitude calculated
according to this formula does not match fig. (5.3) in [1]. If we use
parameters stated in the article proper formula which corresponds to this
figure is

$K_{ij}=\frac{(s-s_{0})}{4m^{2}_{K}}\left[\sum_{p}\frac{f_{i}^{p}f_{j}^{p}}{(s_{p}-s)(s_{p}-s_{0})}+\sum_{n=0}c_{ij}^{n}\left(\frac{s}{4m^{2}_{K}}-1\right)^{n}\right]$
(15)

We believe that it is formula (15) that describes the original $K$ solution of
AMP amplitude.

For the $\hat{M}$ matrix in [1] (3.20) there exists the following
parametrization (so called $M$ solution)

$M_{ij}=\frac{a_{ij}}{(s-s_{0})}+\sum_{p}\frac{{f}^{p}_{i}{f}^{p}_{j}}{(s_{p}-s)}+\sum_{n=0}{c}_{ij}^{n}\left[\frac{s}{4m^{2}_{K}}-1\right]^{n}$
(16)

Again, we believe that there is a mistake in the sign in [1] here and the
proper formula for $M$ solution is

$M_{ij}=\frac{a_{ij}}{(s-s_{0})}-\sum_{p}\frac{{f}^{p}_{i}{f}^{p}_{j}}{(s_{p}-s)}+\sum_{n=0}{c}_{ij}^{n}\left[\frac{s}{4m^{2}_{K}}-1\right]^{n}$
(17)

Note that for $M$ solution parameters of $T$ matrix tends to zero at
$s\to\infty$ while for $K$ solution they tends to unitary limit. The reason of
this behavior is that outside of its scope polynomial background tends to
infinity. We consider the behavior of $M$ solution more appropriate for our
aims and use it. Parameters of original $M$ solution from [1] are listed in
table 1

$s_{0}$ | $s_{1}$ | $f^{1}_{1}$ | $f^{1}_{2}$ | $a_{11}$ | $a_{12}$ | $a_{22}$
---|---|---|---|---|---|---
-0.0074 | 0.9828 | 0.1968 | -0.0154 | 0.1131 | 0.0150 | -0.3216

$c^{0}_{11}$ | $c^{1}_{11}$ | $c^{2}_{11}$ | $c^{3}_{11}$ | $c^{4}_{11}$
---|---|---|---|---
0.0337 | -0.3185 | -0.0942 | -0.5927 | 0.1957
$c^{0}_{12}$ | $c^{1}_{12}$ | $c^{2}_{12}$ | $c^{3}_{12}$ | $c^{4}_{12}$
-0.2826 | 0.0918 | 0.1669 | -0.2082 | -0.1386
$c^{0}_{22}$ | $c^{1}_{22}$ | $c^{2}_{22}$ | $c^{3}_{22}$ | $c^{4}_{22}$
0.3010 | -0.5140 | 0.1176 | 0.5204 | -0.3977

Table 1: Table of coefficients of original AMP amplitude

$s_{0}$ | $s_{1}$ | $f^{1}_{1}$ | $f^{1}_{2}$ | $a_{11}$ | $a_{12}$ | $a_{22}$
---|---|---|---|---|---|---
-0.0074 | 0.9828 | 0 | 0 | 0.1131 | 0 | -0.3216

$c^{0}_{11}$ | $c^{1}_{11}$ | $c^{2}_{11}$ | $c^{3}_{11}$ | $c^{4}_{11}$
---|---|---|---|---
0.0337 | -0.3185 | -0.0942 | -0.5927 | 0
$c^{0}_{12}$ | $c^{1}_{12}$ | $c^{2}_{12}$ | $c^{3}_{12}$ | $c^{4}_{12}$
0 | 0 | 0 | 0 |  0
$c^{0}_{22}$ | $c^{1}_{22}$ | $c^{2}_{22}$ | $c^{3}_{22}$ | $c^{4}_{22}$
0.3010 | -0.5140 | 0.1176 | 0.5204 |  0

Table 2: Table of coefficients for modified amplitude mAMP

##  4 Our modifications

Our aim is to construct the scattering amplitude $\pi\pi\to\pi\pi$ which is
near to the original AMP amplitude, is smooth in the region $M(\pi\pi)\approx
1\,GeV/c^{2}$ and smoothly tends to zero at $M(\pi\pi)>1.5\,GeV/c^{2}$.

We require the proper behavior at $M(\pi\pi)\approx 1\,GeV/c^{2}$ by setting
to zero the connection with $KK$ channel and the connection with $f_{0}(980)$
pole. To suppress an unphysical shoulder in the original amplitude at
$M(\pi\pi)\sim 1.7\,GeV/c^{2}$ we set to zero coefficients at the 4th degree
of the background polynomial.

The parameters of the modified amplitude are given in the table 2. The
original amplitude is shown in the figure 1, the modified one in the figure 2.
The figures show the values, from left to right: Argand plot of the matrix
element $T_{11}$, namely $\pi\pi\to\pi\pi$, its absolute value and Argand
phase.

Figure 1: Original AMP amplitude

Figure 2: Modified AMP amplitude

##  5 Results

The amplitude mAMP is constructed, suitable to describe $\pi\pi\to\pi\pi$
scattering in $S$ wave with $I=0$ in multiparticle PWA. It is not suitable to
describe $KK$ channel.

It is unitary, smooth in the broad range of $M(\pi\pi)$, is near to AMP $M$
solution for $M(\pi\pi)<1.6\ GeV/c^{2}$ except at $M(\pi\pi)\sim
1\,GeV/c^{2}$, smoothly tends to zero at $M(\pi\pi)\to\infty$. It describes
the threshold behavior of $\pi\pi$ scattering and the broad structure at
$M\sim 1400\,MeV/c^{2}$ according to the data known at the time of writing of
[1]. Narrow resonances $f_{0}(980)$, $f_{0}(1500)$ intentionally are not
described, they should be entered into the analysis separately.

The amplitude is built according to formula (17) with coefficients from table
2. This amplitude was used in [5, 6] and the other publications of the VES
group and also in [7]. An early version of the amplitude was used in [8].

The author thanks the members of the VES group for initiating this work,
useful discussions and use of this amplitude in the analysis.

## References

  * [1] K.L. Au, D. Morgan, M.R. Pennington, Meson dynamics beyond the quark model: Study of final-state interactions. Phys. Rev. D35 1633 (1987)
  * [2] J.D. Hansen et al., Formalism and assumptions involved in partial wave analysis of three-meson systems. Nucl. Phys. B81 403 (1974).
  * [3] S.U. Chung and E. Klempt, A Primer on K -matrix Formalism. BNL preprint Uni-Mainz IP-92-03, (1995)
  * [4] Chung, S.U., Brose, J., Hackmann, R., Klempt, E., Spanier, S. and Strassburger, C. (1995), Partial wave analysis in K-matrix formalism. Ann. Phys., 507, pp 404-430. doi:10.1002/andp.19955070504
  * [5] D.I. Amelin et al., Study of resonance production in diffractive reaction $\pi^{-}A\to\pi^{+}\pi^{-}\pi^{-}$, Phys. Lett. B 356 (1995), 595.
  * [6] I.A. Kachaev, D.I. Ryabchikov, High spin resonances produced in $\pi^{+}\pi^{-}\pi^{-}$ and $\pi^{-}\pi^{0}\pi^{0}$ systems at VES setup, MESON 2018, EPJ Web Conf. 199 (2019) 02025
  * [7] C. Adolph et al., Resonance production and $\pi\pi$ S-wave in $\pi^{-}P\to\pi^{-}\pi^{-}\pi^{+}P_{recoil}$ at 190 GeV/c, Phys. Rev. D 95, 032004 (2017)
  * [8] S.U. Chung et al., Exotic and $q\bar{q}$ resonances in the $\pi^{+}\pi^{-}\pi^{-}$ system produced in $\pi^{-}p$ collisions at 18 $GeV/c$, Phys. Rev. D 65, 072001 (2002)