# Utilizing Voronoi Tessellations to Determine Radial Density Profiles

Veronika Dornan Department of Physics and Astronomy
McMaster University
Hamilton, ON, L8S 4M1  William E. Harris Department of Physics and Astronomy
McMaster University
Hamilton, ON, L8S 4M1

###### Abstract

We have developed a novel method of determining 2D radial density profiles for
astronomical systems of discrete objects using Voronoi tessellations. This
Voronoi-based method was tested against the standard annulus-based method on 5
simulated systems of objects, following known Hubble density profiles of
varying parameters and sizes. It was found that the Voronoi-based method
returned radial density fits with lower uncertainties on the fitting
parameters across all 5 systems compared to the annulus-based method. The
Voronoi-based method also consistently returned more accurate estimates of the
total number of objects in each system than the annulus-based method, and this
accuracy increased with increasing system size. Finally, the Voronoi-based
method was applied to two observed globular cluster systems around brightest
cluster galaxies ESO 444-G046 and 2MASX J13272961-3123237 and the results were
compared to previous results for these galaxies obtained with the annulus-
based method. Again, it was found that the Voronoi-based method returned fits
with lower uncertainties on the fitting parameters, and the total number of
globular clusters returned are within errors of the annulus-based method
estimates, however also with lower uncertainties.

galaxies: star cluster – globular clusters: general

##  1 Introduction

Determining the radial density profiles of a variety of astronomical systems
that consist of discrete objects is a frequently encountered problem within
the field. Prominent examples may include stars within a star cluster (King,
1966; Miocchi et al., 2013), systems of globular clusters (GCs) within their
host galaxy (Harris et al., 2015; Dornan & Harris, 2023), or galaxies within a
larger group or cluster (Nagai & Kravtsov, 2005). Thus, in order to get
accurate estimates of the central concentration, profile shape, and total
number of objects within a system it is necessary to have an accurate fit to
the projected radial density profile.

Observationally, the classic method to determine these radial density profiles
is to first determine the 2-dimensional spatial distribution of objects from
photometric observations, then divide this spatial distribution into a series
of concentric circular or elliptical annuli, effectively binning the objects
as a function of radius. For each of these bins the number of objects and the
area of the bin’s annulus are determined in order to obtain a density value.
These density values are then plotted as a function of radius to which a
density profile can then be fit. This method had been applied observationally
to stars in open star clusters (e.g. Alfonso & García-Varela, 2023; Bisht et
al., 2024; Belwal et al., 2024), stars in GCs (e.g. Knödlseder, 2000;
Dalessandro et al., 2015; Simpson et al., 2017), systems of GCs around
galaxies (e.g. Harris et al., 2015, 2016; Harris & Mulholland, 2017; Dornan &
Harris, 2023; Caso et al., 2024; Ines Ennis et al., 2024), galaxies within
galaxy clusters (e.g. Popesso et al., 2007; Newman et al., 2013; Fassbender et
al., 2014; Tozuka et al., 2021), and compact radio sources at the centre of
the Galaxy (Zhao et al., 2022).

This annulus fitting method is useful in its simplicity, but has some
drawbacks. The systems to which it is applied can have thousands of detected
objects within this spatial distribution (Peng et al., 2008, 2011; Harris et
al., 2016; Dornan & Harris, 2023), which when binned, are reduced to only a
few datapoints, resulting in a loss of information about the system, and
uncertainty in the fit. In addition, the number of these binned annuli, and by
extension the radial width they encompass, is arbitrary, which can introduce
further uncertainty in the resulting fit. Finally, as well the annulus method
prescribes a geometrically regular (circular or elliptical) density
distribution, smoothing over substructure such as satellites or local
overdensities.

This work introduces a different general methodology using Voronoi
tessellations for determining radial density profiles. Voronoi tessellation
plots are constructed from a spatial distribution of points, and each
tessellation (Voronoi cell) consists of the area of the distribution that is
closer to the given point than the other points (Okabe et al., 1992). The
result is a distribution of polygons each containing one of the objects from
the input spatial distribution. See Figure 1 for a simple example of a 2D
Voronoi tessellation plot.

The tessellation process allows for construction of the density distribution,
since the area of a Voronoi cell is directly related to the local number
density of objects surrounding the cell. Higher density regions will produce
Voronoi cells with smaller areas, and lower density regions will produce
Voronoi cells with larger areas. Thus, once a Voronoi tessellation plot has
been created from a spatial distribution of objects, the area of each cell can
be inverted to obtain a density value for each object. Dividing the region
into Voronoi cells is a binning process like the standard annular method, but
it essentially carries the definition of bin size to its minimum possible
extreme, with one object per bin. In addition, unlike the classic annulus
method, the Voronoi method makes no assumptions about any geometric regularity
in the spatial distribution.

It should be noted however, that in order to determine surface density using
both annuli and Voronoi tessellations, a distribution of discrete objects is
required. Neither of these methods can be applied to, for example, surface
brightness profiles.

Voronoi tessellations have been applied to the determination of spatial
densities in astronomical research in the past, as they can be very useful in
describing the spatial distribution of observed galaxies, which can exist in
areas of over and under densities (van de Weygaert & Icke, 1989; Icke & van de
Weygaert, 1991). Voronoi tessellations have also been used in algorithms to
automatically detect these over-densities, galaxy clusters, from 2D galaxy
field spatial distributions (Ramella et al., 2001).

In the present paper, we test this Voronoi-cell method against simulated
radial distributions of populations following a simple input distribution (a
Hubble model profile), and compare the results with the standard annular-bin
method. The data used to test this new Voronoi method will be discussed in
section 2. This section will describe the simulated 2D Hubble distributions
and test the method’s ability to return the input radial density profile,
compared to the annulus method. Section 3 will further describe the Voronoi
method in detail and the statistical tests used to compare it to the annulus
method. Section 4 will summarize the results of these statistical tests and
outline the benefits of the Voronoi method for determining radial density
profiles. Finally, section 5 will provide a sample application to observed
systems and discuss future applications of this method and improvements that
can be made.

##  2 Data

In order to test the ability of both the annulus fitting method and the
Voronoi fitting method to accurately return the true radial density profile of
a 2D spatial distribution of objects, a series of artificial systems were
created with known Hubble profiles, defined by

$\sigma=\sigma_{0}\Big{(}1+\frac{r}{r_{0}}\Big{)}^{-b}\,$ (1)

where $\sigma_{0}$ is the density at $r=0$, $r_{0}$ is the core radius, and
the exponent $b$ sets the power-law slope of the profile at large radius.

We generated simulated systems with five different total populations, The
properties of which are detailed in Table 1 (The scale units are arbitrary,
but are simply labelled as ‘pixels’ for convenience). The spatial distribution
of objects for each of these systems was generated with 30 different random
seeds each, allowing for each system to be tested with the same underlying
Hubble distribution but varying positions of objects.

The parameters for the artificial systems were chosen building from our
previous work on the estimated properties of a sample of observed GC systems
hosted by massive elliptical galaxies. These galaxies have already had the
radial density profiles of their GC systems determined with the annulus
fitting method in Dornan & Harris (2023), hereafter referred to as Paper I.
The parameters of the simulated systems were chosen to reflect the range of
profiles seen in a real-life application of this methodology.

It should be noted that in this initial study the Voronoi method is tested
only on circularly symmetric systems, although does not make any assumptions
on the underlying symmetry of the system. This decision was made because the
annulus method, to which we are directly comparing the Voronoi method, does.
In followup work, the Voronoi method will be extended to systems where one or
more satellites are present, with their own subpopulations of objects, in
order to investigate the technique’s ability to handle asymmetrical systems.

Table 1: Artificial System Hubble Profile Parameters

System Number | Objects | Core (pixels) | $b$ |  $r_{max}$ (pixels) |  $\sigma_{o}$ (pixels-2)
---|---|---|---|---|---
(1) | (2) | (3) | (4) | (5) | (6)
1 | 300 | 75 | 1.5 | 3250 | $8.83\times 10^{-4}$
2 | 2000 | 150 | 0.96 | 3250 | $6.93\times 10^{-4}$
3 | 4500 | 150 | 0.96 | 3250 | $1.56\times 10^{-3}$
4 | 7000 | 250 | 1.09 | 3250 | $2.05\times 10^{-3}$
5 | 20000 | 250 | 1.09 | 3250 | $5.86\times 10^{-3}$

Key to columns: (1) Artificial system identification; (2) number of objects in
the system; (3) size of the core of the Hubble profile in units of pixels; (4)
exponent of the Hubble profile; (5) maximum radius objects can be from the
centre of the spatial distribution in units of pixels; (6) density at $r=0$ of
the Hubble profile in units of objects per pixel2.

##  3 Methods

Figure 1:  Left: The 2D spatial distribution of objects with concentric annuli
drawn out from the centre of the system. These annuli are used to bin the
objects and determine densities as a function of radius. Right: The Voronoi
Tessellation plot for the same 2D distribution of objects, with one object per
cell. For both plots cell colour represents comparative density on the same
scale.

###  3.1 Annulus Fitting Method

As noted above, annulus fitting is a standard binning method used to determine
radial density profiles of a variety of circularly and spherically distributed
systems, both simulated and observed (King, 1966; Nagai & Kravtsov, 2005;
Miocchi et al., 2013; Dornan & Harris, 2023; Alfonso & García-Varela, 2023).

An example of a 2D spatial distribution of objects divided into concentric
annuli can be seen in Figure 1, which is for simulated system number 3. These
binned densities can then be plotted as a function of radius and a simple
fitting function such as a power law or Sersic profile can be fit to it.

For these simulated systems the number of annuli was always set to 26 and the
bin widths determined from the range of radii encompassing the innermost and
outermost objects. However, often the 1-2 outermost bins had to be removed for
incompleteness, as often the final few annuli would only host a very low
number of objects and areas. This can be seen in Figure 1 in the corners of
the annulus distribution.

###  3.2 Voronoi Fitting Method

The Voronoi fitting method is capable of retaining far more information about
the system than the annulus fitting method through the use of the Voronoi
tessellations. When a Voronoi cell is generated around each object in the
system, the number of datapoints is equal to the number of objects within the
system. A density value for each datapoint, at the radius of each object, is
determined by simply inverting the area of the Voronoi cell.

In practice, the resulting plot of density versus radius for one object per
cell displays the desired profile, but is also dominated by large stochastic
scatter. A useful approach to reduce the scatter is simply to combine adjacent
cells into larger combined ones. We have experimented with grouping 3, 5, or 7
cells together compared to no grouping at all, as shown in Figure 3. Binning
the Voronoi cells to have five objects per cell proved to be an excellent
compromise between significantly reducing scatter and keeping a large number
of independent datapoints.

An alternative approach to grouping cells would be to combine a given cell
with all its adjacent cells (i.e. those whose sides are in direct contact).
However, for a typical Voronoi tesselation, approximately $89\%$ of all cells
have 6 or more adjacent cells (see Table 1 of Suârez-Plasencia et al., 2021).
Thus the majority of cells in the distribution would be binned in groups of 7
or more, which, as previously mentioned, was found to result in too strong an
averaging and greater culling of datapoints.

Combining cells was done by sorting the cells by density, from highest to
lowest. The highest density cell in the list then has the four spatially
closest objects to it identified, and all five are added into a single bigger
cell and removed from the initial list. This new cell is assigned an (x, y)
position in the grid that corresponds to the averaged positions of the five
member objects, and an area equal to the combined areas of the original cells.
The rest of the binning continues in the same way for the next highest density
individual cell in the list. This procedure ensures that the objects that are
being binned together are spatially close to one another and prevents
situations where the last two objects are on opposite sides of the system, or
oddly shaped (“gerrymandered”) combined cells are created.

Figure 2:  A histogram of the normalized density of Voronoi cells for a Hubble
distribution. The red curve is the log-normal fit to the histogram and the red
line is the average for the Hubble distribution. The blue curve is the fit to
the histogram for a uniform distribution and the blue line is average for the
uniform distribution. The black line is the mean normalized number density,
where $x=n/\langle n\rangle=1$ when $n=\langle n\rangle$ for all distributions

.

Binning the cells together not only reduces scatter, but also reduces a
particular statistical bias that is associated with Voronoi tessellations. A
tessellation of a random, _spatially uniform_ 2D distribution of objects can
be described (Jarai-Szabo & Neda, 2004) as

$f_{2D}(y)=\frac{343}{15}\sqrt{\frac{7}{2\pi}}y^{5/2}\exp{\Big{(}-\frac{7y}{2}\Big{)}}$
(2)

where $y$ is the normalized cell area, $S$, as $y=S/\langle S\rangle$. This
distribution for cell size can be transformed into a distribution for a number
density ($n$) distribution for each cell. Here we will take $x=n/\langle
n\rangle=1/y$ to be the normalized number density. The 2D number density
distribution for a uniform spatial distribution of objects is

$f_{2D}(x)=\frac{343}{15}\sqrt{\frac{7}{2\pi}}x^{9/2}\exp{\Big{(}-\frac{7}{2x}\Big{)}}$
(3)

Figure 3:  Sample radial profile showing the radial density distributions for
1, 3, 5, and 7 objects per cell (going from top down). All annulus data and
fitting is identical, as it is independent from the Voronoi fitting.

For this distribution the mean value of the probability density function does
not equal the mean normalized number density, $x=n/\langle n\rangle$. It is
instead 1.4 times the mean normalized number density, indicating that, for a
uniform spatial distribution, Voronoi tessellations tend to overestimate the
influence of high density areas for the overall system. This can be seen in
Figure 2 when comparing the blue curve, which corresponds to Equation 3, and
the blue vertical line, which represents the curve’s mean, to the black line
which is the expected mean value, where $n=\langle n\rangle$ and $x=1$.

However, for the purposes of this work we are not concerned with the number
density bias for a uniform distribution, but rather for a distribution
following a power-law-like radial profile. When the normalized number density
distribution is plotted as a histogram for the artificial systems in this
work, the distribution can instead be described by a log-normal distribution,
taking the form:

$f_{2D}(x)=\frac{1}{xa\sqrt{2\pi}}\exp{\Big{(}\frac{(\log{x}-b)^{2}}{2a^{2}}\Big{)}}$
(4)

Where $a$ and $b$ are free parameters. This can again be seen in Figure 2,
where the red curve corresponds to Equation 4 and the red vertical line
represents this curve’s mean.

For a Hubble distribution, the average number density of the cells, after
binning into groups of five, is actually 0.90 times the mean normalized number
density. This statistical underestimation of the average number density of the
Voronoi tessellations was found to be small enough that it has a negligible
effect on the radial density profile fits. As will be seen below, the over-
and under-estimations of the density profile at various radii due to simple
randomness across the systems studied were consistently greater than the
influence of this statistical bias, so at this time it is ignored.

##  4 Results

Figure 4:  The density data for the same system using the annulus fitting
method (red datapoints) and the Voronoi fitting method (blue datapoints). Both
datasets are fit with a powerlaw function and compared to the “true” Hubble
profile (green line). Table 2: Bootstrapped Powerlaw Fits with Standard Errors

System Number | Annulus Method | Voronoi Method
---|---|---
1 | $\sigma=(1.131\pm 0.110)r^{(-1.360\pm 0.021)}$ | $\sigma=(1.184\pm 0.084)r^{(-1.370\pm 0.020)}$
2 | $\sigma=(0.852\pm 0.035)r^{(-0.794\pm 0.010)}$ | $\sigma=(0.849\pm 0.029)r^{(-0.784\pm 0.007)}$
3 | $\sigma=(1.899\pm 0.040)r^{(-0.782\pm 0.006)}$ | $\sigma=(1.704\pm 0.012)r^{(-0.753\pm 0.002)}$
4 | $\sigma=(2.637\pm 0.127)r^{(-0.757\pm 0.010)}$ | $\sigma=(2.765\pm 0.011)r^{(-0.758\pm 0.001)}$
5 | $\sigma=(7.377\pm 0.244)r^{(-0.747\pm 0.008)}$ | $\sigma=(8.726\pm 0.013)r^{(-0.780\pm 0.0004)}$

These fitting parameters and standard errors were obtained from 5000
iterations of bootstrapping for both methods.

For each run of each system the density profiles were determined both with the
annulus method and the Voronoi method. This created two different radial
density datasets from the same underlying Hubble distribution. For both the
annular and Voronoi methods, objects within the core of the system (defined in
table 1) were removed. This is because in real observations of systems of
globular clusters the background light intensity within the core of the
central galaxy is typically high enough that identification and measurement of
objects in the core cannot be done (Harris et al., 2015, 2016; Harris, 2023;
Dornan & Harris, 2023; Ines Ennis et al., 2024). Thus, these objects have been
removed for the analysis in order to better re-create the data available for a
real astronomical system, to which this method will be applied in section 5.1.

Figure 1 shows the spatial distribution of the simulated objects and a
visualization of the density fitting method used. Once the objects have been
binned according to the methods used their densities can be plotted as a
function of radius and have profiles fit to them. Here the radii and densities
have been converted from pixels to “arcseconds” according to HST ACS camera
pixel conversions, defined as 1 pixel width $=0.05$ arcseconds (Ryon, 2019),
again to better reflect observed data. The uncertainties associated with the
annulus fitting method density values assume $\sqrt{N}$ statistics.

The uncertainties of the density values from the Voronoi fitting method were
determined from the root mean square scatter around a powerlaw fit to the
densities as a function of radius, for 10 radial bins. This allows the
decreasing uncertainty in cell density to be reflected in the voronoi cells,
as is seen for the annulus bins. Once the uncertainty for each radial bin was
determined it was applied to all density values in that bin.

Once the uncertainties for both methods are determined the Voronoi fitting
method density values were culled to remove outliers. This is done by
determining the average Voronoi density value for the same number of radial
bins as was used for the annuli (although the average Voronoi densities in
these bins are slightly different than the annuli density values). All Voronoi
density values that are within 1.5 standard deviations of these radially
binned means are kept, and all others are culled. The percentage of cells kept
after culling varied for system and run, but was roughly 85%.

In the case of this work, although the distribution was created according to a
Hubble profile, the densities are fit with a simpler powerlaw function due to
the lack of information within the core. These fits can then be compared to
the “true” Hubble profile to determine which method returns the more nearly
correct density profile. An example is shown in Figure 4. If this Voronoi
method was applied to, for example, a single GC, a Hubble profile could then
be used to fit the data as the information in the core would be available.

In order to compare the effectiveness of these two methods of determining
radial profiles there are two important factors to examine: the methods’
precision, and their accuracy. For the precision we compare the uncertainties
of the fits to the data, while for the accuracy, we compare the total numbers
of objects integrated over all radii relative to the input number.

The uncertainties on the fits to the profiles of all systems with both methods
were determined via bootstrapping. The radial density data for both methods
was bootstrapped with 5000 iterations for each system, and powerlaw profile
fitting values were determined, as well as their standard deviations. The
uncertainty for each of these fitting values, the standard error, defined here
in Equation 5, was determined from the standard deviations.

$SE=\frac{\sigma}{\sqrt{n}}$ (5)

Here $\sigma$ is the standard deviation of the scatter of points around the
fitted curve, and $n$ is the length of the dataset used to determine the
fitting parameters.

Figure 5:  A boxplot depicting the results of the number of objects in each
system returned when the radial density profile obtained from each method is
integrated, done across all 30 random object distributions for each system.
The dashed line shows where 100% of all objects are returned. The boxes
represent the interquartile range (IQR), the black lines within the boxes
denote the medians, the diamonds denote the means, and the errorbars encompass
all non-outlier data that lie outside the IQR. The coloured circles represent
outliers, defined as data that lie greater than 1.5 times the IQR past the
upper and lower quartiles

.

In an ideal case, even when the total number of objects in the system is not
completely observed, as in the case of the inner regions of these simulated
systems seen in Figure 1, if the fit is accurate and integrated from zero to
the system’s maximum radius it should return the true number of objects. This
was done for each run of each simulated system, and the results are displayed
in Figure 5 as the total number of objects returned by the fit, relative to
the total input number over the same radial region.

##  5 Discussion

It can be seen from Figure 5 that the Voronoi method is more successful at
returning a total population that is close to the input total. For systems 2
through 5 the Voronoi method’s upper quartiles enclose or very nearly enclose
the true number of objects. For these same systems however, the lower
quartiles of the annulus method lie much further from the true number of
objects, and for systems 3 through 5 the true number of objects is not even
within their lower range. This shows that, certainly for larger systems of
2000 objects or more, the Voronoi method will more consistently and more
accurately return the true number of objects.

Across Figure 5 we can see that the spread in both methods decreases with
increasing system size. Thus, for system 1, the smallest system in this work,
the spread for both methods is significantly higher than for the other
systems. But even for the smallest system the Voronoi method still returns a
slightly more accurate result, with both its mean and median lying closer to
the true number of objects returned than for the annulus method (0.949 and
0.947 compared to 1.033 and 1.047, respectively).

The radial profile fits for each system using each method are listed in Table
2. The precision of each of these methods is compared via the uncertainties on
the fitting parameters for each method. As can be seen in Table 2, the Voronoi
method consistently returns lower uncertainties on both powerlaw fitting
parameters for all systems studied. As system size increases so too does the
difference in method precision, with the parameter uncertainties for system 1
being comparable between methods, to systems 4 and 5 having uncertainties
differing by an order of magnitude.

Figure 6:  The density data for the observed GCSs MASX J13272961-3123237 (top)
and ESO 509-G067 (bottom) using the annulus fitting method (red datapoints)
and the Voronoi fitting method (blue datapoints). Both datasets are fit with a
powerlaw function. The gray shaded region is unobservable due to high surface
brightness.

###  5.1 An Application to Observed Distributions

Here we present a comparison of these two methods when applied to real,
observed GCSs, for which a definition of the “true” underlying radial density
profile is unknown. A common problem in many observed cases is that giant
galaxies may have nearby satellites with their own GCSs, which can overlap the
GCS of the central giant (Dornan & Harris, 2023; Ines Ennis et al., 2024; Lim
et al., 2024). We will return to this separate problem in followup work. For
the purposes of this paper we will focus on two observed examples with minimal
satellite interference. The GCSs chosen for this comparison are those hosted
by the BCGs 2MASSX J13272961-3123237 and ESO 509-G067. Their GCS radial
density profiles have already been determined with the annulus method in Paper
I, and have few significant satellite galaxies. The comparison of these two
methods, applied in the identical way as for the simulated systems, can be
seen in Figure 6. From the work done in Paper I 2MASSX J13272961-3123237 is
estimated to have a GCS population of 9100, and ESO 509-G067 to have a GCS
population 4200, putting them in terms of size close to systems 3 and 4
simulated in this paper.

The radial density fits to these observed GCSs, and the uncertainties on these
fits, using both methods are shown in table 3. The number of GCs in the GCS
and the total mass of the GCS obtained from using these fits are shown in
table 3 as well, using the same methods and photometric completion corrections
as Paper I.

While no statement can be said about the comparative accuracies of these two
methods when applied to observed systems, it can be seen in Figure 6 that the
fits agree well. This is as expected as for ESO 509-G067, as for a system of
this size it was found in Table 2 that the fitting parameters between both
methods agree within uncertainties. However, the fits between methods for
2MASSX J13272961-3123237 agree strikingly well, and is not something that is
expected from the simulations for a system of its size. The further
application of these methods to extremely massive observed systems is needed,
and plans for this is discussed in Section 6.1.

The estimates of the populations themselves are roughly consistent with what
we would expect based on the results from Figure 5, within their
uncertainties. For the larger system, 2MASSX J13272961-3123237, the annulus
method returned a higher estimate than the Voronoi method, albiet with much
less a dramatic difference than expected, as was seen for system 4. For ESO
509-G067, its size would make it roughly comparable to systems 2 and 3, and
although the Voronoi method actually returned a higher estimate than the
annulus method, within the uncertainties of the two methods, this is still be
consistent with the results from Figure 5.

Overall, for both of the returned GCS values for both galaxies in Table 3 both
methods return estimates within their mutual uncertainties. However, the
Voronoi method still returns estimates with lower uncertainties, indicating
that this method is capable of more precise estimates.

##  6 Conclusions

This work develops a new methodology to determine radial density profiles of a
variety of astronomical systems of discrete objects using Voronoi
tessellations. Our main conclusions are as follows:

  * •

This Voronoi method was tested and compared to the classic annulus method for
five simulated systems of varying sizes and steepness of known Hubble
profiles, each run 30 times with different random seeds. These systems were
chosen for their similarities to observed GCSs.

  * •

It was found that for all five systems the standard errors on the fits found
with the Voronoi method were lower than for the fits found with the annulus
method, indicating that the Voronoi method returns more precise fits
regardless of system.

  * •

It was found that the Voronoi method returned more accurate estimates of the
total population of all five systems than the annulus method. The spread in
potential population estimates decreased for both methods with increasing
system size. For systems of roughly more than 1000 objects the annulus method
overestimated the total population by approximately 10%, while the Voronoi
method consistently returned estimates at or within 3% of the true value.

  * •

When these methods were applied to two observed GCSs it was found that the
radial density profile fits produced by the two methods were very similar and
the total $N_{GC}$ estimates were within both methods’ uncertainties. However,
the uncertainties on both the fits and the $N_{GC}$ estimates for both
observed systems were lower for the Voronoi method than for the annulus
method.

Table 3: Comparison of Annulus and Voronoi Methods on Observed Galaxies Galaxy | Value | Annulus Method | Voronoi Method
---|---|---|---
J13272961 | Powerlaw | $\sigma=(16.1\pm 1.5)r^{(-0.92\pm 0.03)}$ | $\sigma=(16.6\pm 0.8)r^{(-0.93\pm 0.01)}$
$N_{GC}$ | $9055\pm 1173$ | $8885\pm 793$
ESO 509-G067 | Powerlaw | $\sigma=(4.5\pm 0.5)r^{(-0.92\pm 0.03)}$ | $\sigma=(5.2\pm 0.02)r^{(-0.93\pm 0.01)}$
$N_{GC}$ | $4174\pm 726$ | $4632\pm 471$

###  6.1 Future Work

Future work with this new method of radial density profile fitting will
include its application to more observed systems of varying sizes. Currently
the authors of this paper plan to apply this method to the remaining GCSs
studied in Paper I, as well as 17 additional BCGs imaged with the HST.

Further improvements to this method also include the ability to accurately
remove objects associated with nearby systems. As briefly mentioned in section
5.1, observed GCSs often have nearby satellite galaxies with their own GCSs
which can contaminate radial density profile estimates of the target system.
The current method of removal of these contaminating systems is to simply mask
out any objects suspected to not belong to the target system (Dornan & Harris,
2023; Lim et al., 2024), or even to simply leave the contaminating system if
it is small enough to not have a drastic affect on the target system (Harris
et al., 2015; Harris & Mulholland, 2017). This method can lead to objects
associated with the target system being removed, or missing various objects
associated with the contaminating system.

With the increase in spatial density data associated with the Voronoi method
compared to the annulus method, it is possible that the radial density profile
of the contaminating system, and its influence on the estimated profile of the
target system, can be determined and subtracted from the target system’s fit.
This would not only give a far more accurate estimate of the profile for
systems in dense, clustered environments, but would also allow for an estimate
of the smaller contaminating system’s radial density profile. This can be
particularly helpful in situations where the GCS distributions of two
interacting galaxies is trying to be studied (eg. Ines Ennis et al., 2024).
Work on developing this method is in progress.

##  7 Acknowledgements

This work was supported by a Discovery Grant to WEH from the Natural Sciences
and Engineering Research Council of Canada (NSERC).

## References

  * Alfonso & García-Varela (2023) Alfonso, J., & García-Varela, A. 2023, A&A, 677, A163, doi: 10.1051/0004-6361/202346569
  * Belwal et al. (2024) Belwal, K., Bisht, D., Singh Bisht, M., et al. 2024, arXiv e-prints, arXiv:2403.04532, doi: 10.48550/arXiv.2403.04532
  * Bisht et al. (2024) Bisht, D., Bisht, D., Raj, A., et al. 2024, New A, 109, 102205, doi: 10.1016/j.newast.2024.102205
  * Caso et al. (2024) Caso, J. P., Ennis, A. I., & De Bórtoli, B. J. 2024, MNRAS, 527, 6993, doi: 10.1093/mnras/stad3602
  * Dalessandro et al. (2015) Dalessandro, E., Ferraro, F. R., Massari, D., et al. 2015, ApJ, 810, 40, doi: 10.1088/0004-637X/810/1/40
  * Dornan & Harris (2023) Dornan, V., & Harris, W. E. 2023, ApJ, 950, 179, doi: 10.3847/1538-4357/accbc3
  * Fassbender et al. (2014) Fassbender, R., Nastasi, A., Santos, J. S., et al. 2014, A&A, 568, A5, doi: 10.1051/0004-6361/201423941
  * Harris (2023) Harris, W. E. 2023, ApJS, 265, 9, doi: 10.3847/1538-4365/acab5c
  * Harris et al. (2016) Harris, W. E., Blakeslee, J. P., Whitmore, B. C., et al. 2016, Astrophysical Journal, 817, 58, doi: 10.3847/0004-637X/817/1/58
  * Harris et al. (2015) Harris, W. E., Harris, G. L., & Hudson, M. J. 2015, ApJ, 806, 36, doi: 10.1088/0004-637X/806/1/36
  * Harris & Mulholland (2017) Harris, W. E., & Mulholland, C. J. 2017, ApJ, 839, 102, doi: 10.3847/1538-4357/aa6a59
  * Icke & van de Weygaert (1991) Icke, V., & van de Weygaert, R. 1991, QJRAS, 32, 85
  * Ines Ennis et al. (2024) Ines Ennis, A., Caso, J. P., & Bassino, L. P. 2024, arXiv e-prints, arXiv:2403.00065, doi: 10.48550/arXiv.2403.00065
  * Jarai-Szabo & Neda (2004) Jarai-Szabo, F., & Neda, Z. 2004, arXiv e-prints, cond, doi: 10.48550/arXiv.cond-mat/0406116
  * King (1966) King, I. R. 1966, AJ, 71, 276, doi: 10.1086/109918
  * Knödlseder (2000) Knödlseder, J. 2000, A&A, 360, 539, doi: 10.48550/arXiv.astro-ph/0007442
  * Lim et al. (2024) Lim, S., Peng, E. W., Côté, P., et al. 2024, arXiv e-prints, arXiv:2403.09926, doi: 10.48550/arXiv.2403.09926
  * Miocchi et al. (2013) Miocchi, P., Lanzoni, B., Ferraro, F. R., et al. 2013, ApJ, 774, 151, doi: 10.1088/0004-637X/774/2/151
  * Nagai & Kravtsov (2005) Nagai, D., & Kravtsov, A. V. 2005, ApJ, 618, 557, doi: 10.1086/426016
  * Newman et al. (2013) Newman, A. B., Treu, T., Ellis, R. S., et al. 2013, ApJ, 765, 24, doi: 10.1088/0004-637X/765/1/24
  * Okabe et al. (1992) Okabe, A., Boots, B., & Sugihara, K. 1992, Spatial tessellations. Concepts and Applications of Voronoi diagrams (Wiley)
  * Peng et al. (2008) Peng, E. W., Jordán, A., Côté, P., et al. 2008, ApJ, 681, 197, doi: 10.1086/587951
  * Peng et al. (2011) Peng, E. W., Ferguson, H. C., Goudfrooij, P., et al. 2011, ApJ, 730, 23, doi: 10.1088/0004-637X/730/1/23
  * Popesso et al. (2007) Popesso, P., Biviano, A., Böhringer, H., & Romaniello, M. 2007, A&A, 464, 451, doi: 10.1051/0004-6361:20054708
  * Ramella et al. (2001) Ramella, M., Boschin, W., Fadda, D., & Nonino, M. 2001, A&A, 368, 776, doi: 10.1051/0004-6361:20010071
  * Ryon (2019) Ryon, J. E. 2019, Advanced Camera for Surveys Instrument Handbook for Cycle 27 v. 18.0
  * Simpson et al. (2017) Simpson, J. D., De Silva, G., Martell, S. L., Navin, C. A., & Zucker, D. B. 2017, MNRAS, 472, 2856, doi: 10.1093/mnras/stx2174
  * Suârez-Plasencia et al. (2021) Suârez-Plasencia, L., Herrera-Macîas, J. A., Legôn-Pêrez, C. M., et al. 2021, in Lecture Notes of the Institute for Computer Sciences, Social Informatics and Telecommunications Engineering, doi: 10.1007/978-3-030-69839-3_13
  * Tozuka et al. (2021) Tozuka, M., Kawabata, K., Fukazawa, Y., Nishiura, S., & Hiraga, J. 2021, Stars and Galaxies, 4, 6, doi: 10.32231/starsandgalaxies.4.0_6
  * van de Weygaert & Icke (1989) van de Weygaert, R., & Icke, V. 1989, A&A, 213, 1
  * Zhao et al. (2022) Zhao, J.-H., Morris, M. R., & Goss, W. M. 2022, ApJ, 927, L6, doi: 10.3847/2041-8213/ac54be


# How prominent are microemulsion phases in 2D electron systems?

Sandeep Joy  Brian Skinner  Department of Physics, Ohio State University,
Columbus, OH 43210, USA

###### Abstract

In two-dimensional electronic systems, direct first-order phase transitions
are prohibited as a consequence of the long-range Coulomb interaction, which
implies a stiff energetic penalty for macroscopic phase separation. A
prominent proposal is that any direct first-order transition is instead
replaced by a sequence of “microemulsion” phases, in which the two phases are
mixed in patterns of mesoscopic domains. In this note, we comment on the range
$\Delta n$ of average electron density that such microemulsion phases may
occupy. We point out that, even without knowing the value of a
phenomenological parameter associated with surface tension between the two
phases, one can place a fairly strong upper bound on the value of $\Delta n$.
We make numerical estimates for $\Delta n$ in the case of the Fermi liquid to
Wigner crystal transition and find $\Delta n$ to be on the order of $10^{7}$
cm-2. This value is much smaller than the width of the phase transition
observed in experiments, suggesting that disorder is a more likely explanation
for the apparent broadening of the transition.

Over two decades ago, it was established that there can be no direct first-
order transition between solid and liquid phases of a two-dimensional electron
system (2DES) [1, 2, 3]. The authors of Refs. [1, 2, 3] proposed that the
nominal first-order transition is instead replaced by a series of
“microemulsion” phases, in which the two distinct states are blended together
in spatially-structured mesoscopic domains (similar ideas appeared earlier in
Ref. [4]). Subsequent work has developed this idea further and explored its
experimental implications [5, 6, 7]. The possibility of microemulsion phases
is now commonly invoked when experiments show a phase transition in a 2DES
that occupies a finite window of average electron density rather than an
abrupt transition at a single density (e.g., Refs. [8, 9, 10]). In this short
note, we comment on how large, exactly, can be the window of average density,
$\Delta n$, occupied by microemulsion phases. After deriving a general bound
on $\Delta n$ we focus specifically on the quantum melting transition between
the Wigner crystal and Fermi liquid states.

Figure 1: Schematic depiction of the conventional Maxwell construction for a
first-order phase transition. The red and blue curves represent the free
energy per unit area of phases $1$ and $2$, respectively, as a function of the
average density $\bar{n}$. The common tangent represents the coexistence line,
in which there a (macroscopic) phase separation. During this phase separation
regime, both phases have the same chemical potential (the slope of the
tangent).

A first-order phase transition and its associated (macroscopic) phase
coexistence is typically understood using the framework of the Maxwell
construction. Before considering the role of long-range Coulomb interactions
between phase-segregated domains, we first briefly summarize the conventional
Maxwell construction. Consider a (first-order) phase transition between two
different phases, which we label as 1 and 2. If $\varepsilon_{i}(n)$ denotes
the energy per particle as a function of the particle density $n$ in phase
$i=1,2$, then at zero temperature, the free energy per unit area in phase $i$
is $f_{i}(n)=n\varepsilon_{i}(n)$. In the usual Maxwell construction for phase
coexistence, one considers that the total free energy for a fixed number of
particles is minimized with respect to the area fraction $x$ that is occupied
by one of the two phases (say, phase 2). That is, if $n_{1},n_{2}$ denote the
particle densities in spatial regions occupied by phases $1$ and $2$,
respectively, then the region of average density which produces phase
coexistence is found by minimizing the total free energy
$(1-x)f_{1}(n_{1})+xf_{2}(n_{2})$ with respect to $x$ and $\delta
n=n_{2}-n_{1}$ for a fixed average density $\bar{n}=(1-x)n_{1}+xn_{2}$. From
this procedure, one can easily show that, within the regime of phase
coexistence, the two phases have the same chemical potential
$\mu_{i}=df_{i}(n)/dn$ (see Fig. 1). 111The slope of the common tangent is
given by
$(f_{2}(n_{2})-f_{1}(n_{1}))/(n_{2}-n_{1})=\mu_{2}(n_{2})=\mu_{1}(n_{1})$.
This condition is equivalent to the more widely recognized version of the
Maxwell construction, known as the equal area rule:
$\int_{n_{1}}^{n_{2}}(\mu(n)-\mu^{\star})\,dn=0$, where $\mu^{\star}$ is the
chemical potential of the system during phase coexistence and $\mu(n)$ denotes
the chemical potential of whichever uniform phase has lower energy at density
$n$.

In electron systems, however, there is an additional Coulombic energy cost
associated with the difference between the local electron density ($n_{1}$ or
$n_{2}$) and the charge density $+e\bar{n}$ of the positive background.
(Throughout this note we assume that the gate electrode containing
neutralizing background charge is much further away from the plane of the 2DES
than the size of microemulsion domains; allowing this gate distance $d$ to be
finite truncates the long-ranged Coulomb interaction and gives a lower bound
on the coexistence window $\Delta n$ that is proportional to $1/d$ [1, 2]).
This Coulombic energy suppresses phase coexistence and generally prohibits the
formation of macroscopic phase domains (since a region of size $L$ and charge
density $\eta$ has a Coulombic self-energy per unit area that grows
extensively with the region size, $\sim\eta^{2}L$) [11]. Instead, one ends up
with the microemulsion phase, where stripes or droplets of the less abundant
phase are interspersed periodically among the more abundant phase.

Argument for the absence of a direct first-order transition. Before discussing
an upper limit on $\Delta n$, we first briefly recapitulate the argument that
there can be no direct first-order phase transitions in a 2DES (see Ref. [6]
for a more complete review). The argument proceeds by considering the system
at the critical density $n=n_{c}$ for which the two phases have equal energy
(the nominal location of the first-order transition) and showing that one can
construct a trial microemulsion state with the same average density but lower
free energy than either pure phase. The existence of such a lower-energy state
implies that the system must not pass directly from a pure phase 1 to a pure
phase 2 as a function of increasing average density, and a direct first-order
transition between the two pure phases is precluded. The argument is simplest
if one assumes that the interaction between electrons is of the form
$V(r)=k/r^{\alpha}$, where $1<\alpha<2$. The limit $\alpha\rightarrow 1$ is
discussed below.

Consider a trial state consisting of alternating stripes of phases 1 and 2,
each with the same width $\ell$ (assumed to be much longer than the inter-
electron spacing) and having uniform electron density $n_{1}=n_{c}-\delta n$
(phase 1) and $n_{2}=n_{c}+\delta n$ (phase 2), so that the global average
density is $n_{c}$. The free energy per area of this phase, relative to either
uniform phase at $n=n_{c}$, can be written

$\delta f=f_{\text{bulk}}+f_{\text{Coulomb}}+f_{\text{surface}}.$ (1)

The first contribution, $f_{\text{bulk}}$, represents the change in free
energy (per unit area) associated with shifting some electrons from phase 1 to
phase 2, which has a lower chemical potential. To leading order in $\delta n$,
$f_{\text{stripe}}$ is given by:

$f_{\text{bulk}}=-\frac{\delta n}{2}(\mu_{1}(n_{c})-\mu_{2}(n_{c})).$ (2)

The second contribution, $f_{\text{Coulomb}}$, corresponds to the
electrostatic energy cost associated with the two stripe regions having an
overall net charge $\pm e\delta n$ per unit area:

$f_{\text{Coulomb}}=A(e\delta n)^{2}V(\ell)\ell^{2},$ (3)

where $A$ is a numerical constant of order unity.

The final contribution, $f_{\text{surface}}$, represents the energy required
to create an interface between the two phases:

$f_{\text{surface}}=\frac{\sigma}{\ell},$ (4)

where $\sigma$ is the interfacial tension.

Keeping $\ell$ fixed, we minimize $\delta f$ with respect to $\delta n$, which
gives

$\delta n=\frac{\mu_{1}(n_{c})-\mu_{2}(n_{c})}{4Ake^{2}\ell^{2-\alpha}}.$ (5)

The corresponding (minimum) free energy change is given by:

$\delta
f=-\frac{\left(\mu_{1}(n_{c})-\mu_{2}(n_{c})\right)^{2}}{16Ake^{2}\ell^{2-\alpha}}+\frac{\sigma}{\ell}.$
(6)

Notice that so long as $1<\alpha<2$, the change in free energy $\delta f$ is
negative whenever $\ell$ is sufficiently long. Thus, we have constructed a
trial microemulsion state that has a lower energy than either uniform state,
and there must not be a direct first-order transition between two uniform
states.

In the limit $\alpha\rightarrow 1$, both terms in Eq. 6 are proportional to
$1/\ell$, and the analysis is inconclusive. However, a more carefully
constructed trial state, in which the electron density is allowed to vary
within each stripe, still yields a negative value of $\delta f$ in the limit
of very long $\ell$ [3].

Coulomb-frustrated phase coexistence. Because of the long-ranged Coulomb
interactions between regions of different charge density, phase coexistence
cannot be understood without simultaneously considering the spatial structure
of domains. This problem has been studied extensively for both 2D and 3D
electron gases (3D in [12, 13, 14], 2D in [15], and both in [16, 17]).

Here we briefly recap the primary results of these studies for a 2DES, and
show how the Maxwell construction for phase coexistence is modified. Following
Ref. [15], we assume uniform density within each domain. We comment on the
effect of relaxing this assumption below.

One can generically write the free energy per unit area of a mixed phase as

$f=(1-x)f_{1}(n_{1})+xf_{2}(n_{2})+e_{m},$ (7)

where $e_{m}$ is the energy density of mixing, which contains both
electrostatic and surface energy terms. The energy $e_{m}$ has been considered
in detail in Ref. [15] for both the droplet and the stripe geometry. For the
case of droplet configurations, $e_{m}$ is well-approximated by 222In stripe
configurations, there is an additional logarithmic factor [15].

$e_{m}=\frac{8}{\sqrt{3}}|n_{2}-n_{1}|\sqrt{\frac{e^{2}\sigma}{\epsilon}}x(1-x),$
(8)

where $\sigma$ is the surface tension. (See Appendix A for a brief derivation
of Eq. 8.) The corresponding optimal droplet size $R_{d}$ is

$R_{d}=\frac{\sqrt{e^{2}\sigma/\epsilon}}{\left(e^{2}/\epsilon\right)|n_{2}-n_{1}|}\frac{\sqrt{3}}{2\sqrt{x\left(1-x\right)}}.$
(9)

(A similar expression was presented in the context of quantum Hall stripe and
bubble phases in Ref. [11].)

Figure 2: Schematic of the modified Maxwell construction accounting for long-
range interactions. The phase coexistence curve (purple line) exhibits a
discontinuity in chemical potential between the two phases, as given by Eq.
11. The corresponding window $\Delta n=n_{2}-n_{1}$ of phase coexistence is
smaller than in the conventional Maxwell construction (Fig. 1).

Thus, both the droplet size and the energy of mixing depend on the parameter
$\beta\equiv\sqrt{e^{2}\sigma/\epsilon}$, which has units of energy. One can
think of $\beta$ as a parameterization of surface tension: large $\beta$
favors large domains so that the total amount of surface area in the sample is
smaller, while small $\beta$ favors small domains. Minimizing the free energy
in Eq. 7 with respect to ${\delta n=n_{2}-n_{1}}$ (as in the Maxwell
construction)333Similarly, minimizing the free energy in Eq. 7 with respect to
$x$ produces a relation between the difference in pressure between the two
phases [15]., we arrive at the following relation:

$\mu_{1}(n_{1})-\mu_{2}(n_{2})=\frac{8}{\sqrt{3}}\beta.$ (10)

That is, while in the conventional Maxwell construction the two phases have
the same chemical potential, the inclusion of the term $e_{m}$ implies that
two coexisting phases have a difference in chemical potential $\sim\beta$.444
We also note that minimizing Eq. 7 gives densities $n_{1}$ and $n_{2}$ that
vary as a function of the area fraction $x$, unlike the conventional Maxwell
construction, in which the densities of the two phases remain fixed as the
average density is varied. This difference in chemical potential can be
interpreted as a result of electrochemical equilibrium, in which the
difference in electrostatic potential between neighboring domains must be
compensated by a difference in the intrinsic chemical potentials.

Upper bound on $\Delta n$. Having outlined the modified Maxwell construction,
we now consider the bound it places on the window of phase coexistence. The
previous analysis provides two fundamental criteria for the presence of a
microemulsion phase over any finite interval $\Delta n$ of the average density
$\bar{n}$. First, there is a maximum value for the difference in chemical
potential between the two phases, effectively determining an upper limit for
$\beta$. In a scenario for which phase coexistence is barely possible (i.e.,
it exists over a very small window of $\bar{n}$), the two phases have
densities that are very close to the value $n_{c}$ at which the two phases
have the same energy, $f_{1}(n_{c})=f_{2}(n_{c})$. In this scenario, the
difference in chemical potential is near its maximum value,
$\Delta\mu=\mu_{1}(n_{c})-\mu_{2}(n_{c})$. Correspondingly phase coexistence
is only permitted when

$\beta<\frac{\sqrt{3}}{8}\Delta\mu.$ (11)

Equation 11 is derived under the assumption that each domain has a spatially
uniform density. In fact, if $\beta$ is large compared to $\Delta\mu$ then
this uniform density approximation breaks down [16, 3], and one should instead
do a more careful calculation that allows for the charge density within each
domain to increase near the domain walls. For large $\beta/\Delta\mu$, this
calculation can produce a stable microemulsion phase, but only over a window
of density that is exponentially small in the parameter $\beta/\Delta\mu$. We
can thus conclude that in order for there to be a significant window of phase
coexistence, $\beta$ must be small enough that Eq. 11 must not be strongly
violated.

On the other hand, if $\beta$ is too small, then the length scale associated
with microemulsion droplets becomes unrealistically small. Specifically, phase
coexistence is not possible unless the corresponding minimal droplet size (the
minimum value of $R_{d}$ is at $x=1/2$ in Eq. 9) is much larger than the
inter-electron spacing. That is, $R_{d}\gg n_{1,2}^{-1/2}\approx
n_{c}^{-1/2}$. This condition produces the strong inequality

$n_{c}^{-1/2}\ll\frac{\sqrt{3}\,\beta}{\left(e^{2}/\epsilon\right)\left|n_{2}-n_{1}\right|}.$
(12)

Equation 12 sets a _lower_ bound on the parameter $\beta$:

$\frac{e^{2}|n_{2}-n_{1}|}{\sqrt{3}\epsilon n_{c}^{1/2}}\ll\beta.$ (13)

Combining Eq. 11 with Eq. 13 produces an upper bound for the difference in
density between the two phases:

$\frac{\Delta
n}{n_{c}}\ll\frac{3}{8}\frac{\Delta\mu}{e^{2}n_{c}^{1/2}\big{/}\epsilon},$
(14)

Equation 14 constitutes the main result of this note.

We emphasize that Eq. 14 makes no assumptions about the value of the surface
tension parameter $\beta$. In a given system, the size $\Delta n$ of the
microemulsion window of course depends on $\beta$, but even without knowing
the value of $\beta$ one can place a strong upper bound on $\Delta n$ by
combining the requirements that the system should reach electrochemical
equilibrium and that droplets should be larger than the interelectron spacing.

Wigner crystal to Fermi liquid transition. We now use the inequality in Eq. 14
to make estimates for the range of phase coexistence in the well-studied case
of the quantum melting of the Wigner crystal (WC) phase. For this estimate we
use the numerically fitted formulas presented in Ref. [18], derived from
quantum Monte Carlo methods, for the energy per electron of the FL and WC
phases. The resulting estimated upper bounds on $\Delta n/n_{c}$ are given in
Table 1. Ref. [18] considers WC trial states with both ferromagnetic (FM) and
antiferromagnetic (AFM) spin ordering, and suggests a narrow window of the
interaction parameter $r_{s}$ for which the AFM WC has the lowest energy. For
completeness we provide estimates of $\Delta n/n_{c}$ for all three possible
transitions. A similar analysis is presented for the case of phase transitions
between bubble phases in the Landau levels of graphene in Appendix B. In that
case we find that microemulsion phases occupy a range of filling factors
$\delta\nu$ of order $10^{-2}$.

Transition | Critical $r_{s}$ | $\max\Delta n/n_{c}$
---|---|---
PM FL to FM WC | $32.89$ | $1.16\times 10^{-3}$
PM FL to AF WC | $31.26$ | $8.45\times 10^{-4}$
AF WC to FM WC | $37.40$ | $3.04\times 10^{-4}$
Table 1: The bounds on $\Delta n/n_{c}$ for various WC-FL transitions are
tabulated using Eq. 14 and numerical results for the chemical potential of
each phase from [18]. Here, FM, AF, and PM represent the Ferromagnetic,
Antiferromagnetic, and Paramagnetic phases, respectively.

In Table II we estimate the value of $\Delta n$, in absolute units, for
different 2DES systems for which a WC-FL transition has been reported (see
Table. 2). In each case, $\Delta n$ is of order $10^{7}$ cm-2, which is
generally too small to be observed in current experiments.

System |  $n_{c}$ (cm-2) | Ref |  $\Delta n$ (cm-2)
---|---|---|---
MgZnO/ZnO | $1.6\times 10^{10}$ | [19] | $1.85\times 10^{7}$
AlAs | $2\times 10^{10}$ | [20] | $2.32\times 10^{7}$
GaAs (holes) | $0.77\times 10^{10}$ | [21] | $0.89\times 10^{7}$
Table 2: The maximum density window of microemulsion phase coexistence between
WC and FL phases is estimated for different experimental 2DES platforms. The
listed value of $\Delta n$ corresponds to the PM FL to FM WC transition (the
first row in Table. 1).

Disorder-induced density fluctuation. So far we have been considering a
nondisordered 2DES, and we have shown that the maximum value of $\Delta n$ for
such systems is very small, on the order of $10^{7}$ cm-2 for the WC-FL
transition. In real experiments, disorder tends to produce a much larger
modulation of the density even in the cleanest samples. The finite width of
the WC-FL transition seen in experiments (e.g., Refs. [19, 20, 21, 22]) is
therefore more likely to be attributable to disorder rather than to
microemulsion physics.

Specifically, disorder creates a spatially-varying electric potential, so that
in equilibrium the local chemical potential also varies spatially. When the
electron density is close to $n_{c}$, the resulting modulation in chemical
potential produces puddles of the minority phase embedded in the majority
phase [23, 24, 25, 26], which appears experimentally as an effective smearing
of the transition. For example, a common source of disorder in experimental
samples is a three-dimensional concentration $N_{\text{imp}}$ of impurity
charges embedded in the substrate, which produces relatively large and long-
wavelength fluctuations of the 2DES density even when $N_{\text{imp}}$ is low
[27]. One can estimate the corresponding value of $\Delta n$ using a self-
consistent a Thomas-Fermi approximation [24], which gives $\Delta
n_{\text{imp}}\sim(e^{2}N_{\text{imp}}\big{/}\epsilon\,|dn/d\mu|)^{1/2}$ [27].
In the case of a 2DES at large $r_{s}$, $|dn/d\mu|$ is of the order of
$\epsilon\,n_{c}^{1/2}/e^{2}$, resulting in $\Delta
n_{\text{imp}}\sim(N_{imp}^{2}n_{c})^{1/4}$. In the case of MgZnO/ZnO, for
example, $\Delta n_{\text{imp}}$ was estimated to be on the order of several
times $10^{9}$ cm-2 even for a bulk impurity concentration as low as $10^{14}$
cm-3 ($\sim 10$ parts per billion) [19].

Closing remarks. We have shown how one can arrive at an upper bound for the
window $\Delta n$ of phase coexistence by combining two criteria: (i) the
existence of electrochemical equilibrium between the two phases, which puts a
lower bound on the surface tension parameter $\beta$; and (ii) the requirement
that the size of a domain be larger than the interelectron spacing, which puts
an upper bound on $\beta$. Combining these two inequalities produces Eq. 14,
which is a seemingly generic upper bound on $\Delta n$ that depends only on
the difference in chemical potential $\Delta\mu$ between the two phases at the
critical point. (An equivalent expression for the context of quantum Hall
phase transitions is presented as Eq. 19 in Appendix B.) For the case of the
WC-FL transition, we use numerical results for $\Delta\mu$ [18] to arrive at
an upper bound $\Delta n/n_{c}\lesssim 10^{-3}$. This very small window of
density is consistent with quantum Monte Carlo studies, which have failed to
find evidence for a regime of phase coexistence at the WC-FL transition [18,
28].

We emphasize that the inequality we derive only applies to phase coexistence
in the sense of mesoscopic domains. That is, one can only talk about “phase 1
coexisting with phase 2” if the domains of each phase are much larger than the
interelectron spacing, so that locally each region resembles the corresponding
bulk phase. However, we cannot rule out the possibility that an instability
toward microscopic level mixing (a violation of Eq. 12) implies the existence
of a new phase in the vicinity of $n_{c}$, distinct from either of the bulk
phases 1 or 2. (For example, one can think of the quantum Hall stripe and
bubble phases [29, 30, 11] as similar to a “microemulsion” phase with a
modulation of density at the scale of the interelectron spacing.) Such new
phases cannot be described by a Maxwell construction type argument.

Acknowledgments: The authors are grateful to Steve Kivelson and Boris Spivak
for helpful discussions and to Andrea Young and Fangyuan Yang for a related
collaboration. This work was supported by the NSF under Grant No. DMR-2045742.

##  Appendix A Electrostatic energy of a microemulsion droplet

Here we provide a calculation of the mixing energy for a microemulsion phase
in which the minority phase (which we take to be phase 2) forms droplets
embedded in the majority phase (phase 1). To model this system, we use a
Wigner-Seitz cell approximation (see Fig. 3), dividing the entire sample into
approximately circular cells of radius $R_{c}$. The total number of cells is
$N_{c}=A/\pi R_{c}^{2}$, where $A$ is the total area of the sample and $R_{d}$
is the radius of the inner disk. The area fraction of phase 2 is represented
by $x=(R_{d}/R_{c})^{2}$. The droplet geometry is applicable when $x$ is close
to either $0$ or $1$. When $x$ is close to $1/2$, a stripe geometry has lower
energy. Nevertheless, one can use the droplet calculation to arrive at an
approximate formula that interpolates across the entire range of $x$ [15], as
we discuss further below.

To ensure that each cell is charge neutral, we introduce a uniformly charged
background with a density of $-\bar{n}$ (where $\bar{n}=(1-x)n_{1}+xn_{2}$).
Within this model, the local density is given by:

$\delta n(\rho)=\left\\{\begin{array}[]{cc}\delta n_{2}&0<\rho<R_{d}\\\ \delta
n_{1}&R_{d}<\rho<R_{c},\end{array}\right.$ (15)

where $\rho$ is the radial coordinate in the plane, $\delta n_{2}\equiv
n_{2}-\bar{n}$, and $\delta n_{1}=n_{1}-\bar{n}$.

We evaluate the electrostatic energy using the electrostatic potential
$V(\rho)$ of a uniform charge disk of radius $R$ and charge density $\sigma$
[31]:

$V(\rho)=\left\\{\begin{array}[]{cc}\frac{4\sigma}{\epsilon}E\left[\left(\frac{\rho}{R}\right)^{2}\right],&0<\rho<R\\\
\frac{4\sigma}{\epsilon}\left(\left(\frac{\rho}{R}\right)E\left[\left(\frac{\rho}{R}\right)^{2}\right]+\frac{\left(1-\left(\rho/R\right)^{2}\right)}{\left(\rho/R\right)}K\left[\left(\frac{\rho}{R}\right)^{2}\right]\right),&R<\rho<\infty\end{array}\right.$
(16)

Figure 3: Diagram representing the construction of Wigner-Seitz cell. The
background density $-\bar{n}<0$ ensures the charge neutrality of the entire
cell. Inside the inner disk, the total density is given by $\delta n_{2}>0$.
In the outer shell, the total density is given by $\delta n_{1}<0$.

The electrostatic energy of a droplet is then given by:

$E_{c}=\frac{8\pi e^{2}(n_{2}-n_{1})^{2}R_{c}^{3}}{3\epsilon}f(x)$ (17)

where $f(x)$ is the function

$f(x)=x\left(\sqrt{x}+x-(1+x)E(x)+(1-x)K(x)\right).$ (18)

Here, $K(x)$ and $E(x)$ are the complete elliptic integrals of the first and
second kind, respectively.

As mentioned above, the configuration in which drops of phase 2 are embedded
in phase 1 appears only when $x$ is sufficiently close to zero. At small $x$,
Eq. 18 gives $f(x)\simeq x^{3/2}$. At $x$ close to $1$ the roles of phase 1
and 2 are reversed, implying $f(x)\simeq(1-x)^{3/2}$. Following the results of
Ref. [15], these two limits imply an approximate expression
$f(x)\approx\left(x(1-x)\right)^{3/2}$ for the entire range of $x$.

By minimizing the sum of the total Coulomb energy ($N_{c}\times E_{c}$) and
the total surface energy ($N_{c}\times 2\pi\sigma R_{d}$) with respect to
$R_{d}$, we obtain the result given in Eq. 8 and E. 9.

##  Appendix B Density range of microemulsion for bubble phases: Illustration
using the $n=3$ Landau level in graphene

Here, as an example of the generality of our result, we use Eq. 14 to place a
bound on the range of filling factors associated with the phase coexistence
between the so-called bubble phase in the $n=3$ Landau Level [29, 30].

Transition | $\nu_{c}$ |  $\max\Delta\nu/\nu_{c}$ (Eq. 19)
---|---|---
$M=1$ to $M=2$ | $0.18$ | $2.25\times 10^{-2}$
$M=2$ to $M=3$ | $0.28$ | $2.56\times 10^{-2}$
Table 3: The bounds on $\Delta\nu/\nu_{c}$ for the phase transition between
bubble phases in the $n=3$ Landau levels of graphene [32]. The integer $M$
labels the number of electrons per bubble.

Inserting units of the magnetic length, we arrive at an alternate formulation
of Eq. 14:

$\Delta\nu\ll\frac{3\sqrt{2\pi}}{8}\frac{\Delta\mu(\nu_{c})}{e^{2}/\epsilon
l_{B}}\nu_{c}^{1/2}.$ (19)

In this equation, $\Delta\nu$ denotes the interval of filling factors
corresponding to phase coexistence, $\nu_{c}$ signifies the partial filling
factor at which the energies of the two phases are equal, and
$\Delta\mu(\nu_{c})$ represents the discrepancy in chemical potential between
the phases at $\nu=\nu_{c}$. As an example, the results for $n=3$ Landau
levels in graphene are presented in Table 3, with results for the energy per
electron taken from Ref. [32].

## References

  * Spivak [2003] B. Spivak, Phase separation in the two-dimensional electron liquid in MOSFET’s, Phys. Rev. B 67, 125205 (2003).
  * Spivak and Kivelson [2004] B. Spivak and S. A. Kivelson, Phases intermediate between a two-dimensional electron liquid and wigner crystal, Phys. Rev. B 70, 155114 (2004).
  * Jamei  _et al._ [2005] R. Jamei, S. Kivelson, and B. Spivak, Universal aspects of coulomb-frustrated phase separation, Phys. Rev. Lett. 94, 056805 (2005).
  * Kivelson and Trugman [1986] S. Kivelson and S. A. Trugman, Quantization of the hall conductance from density quantization alone, Phys. Rev. B 33, 3629 (1986).
  * Spivak and Kivelson [2005] B. Spivak and S. A. Kivelson, Drag resistance of two-dimensional electronic microemulsions, Phys. Rev. B 72, 045355 (2005).
  * Spivak and Kivelson [2006] B. Spivak and S. A. Kivelson, Transport in two dimensional electronic micro-emulsions, Annals of Physics 321, 2071 (2006).
  * Spivak  _et al._ [2010] B. Spivak, S. V. Kravchenko, S. A. Kivelson, and X. P. A. Gao, Colloquium: Transport in strongly correlated two dimensional electron fluids, Rev. Mod. Phys. 82, 1743 (2010).
  * Hatke  _et al._ [2015] A. T. Hatke, Y. Liu, L. W. Engel, M. Shayegan, L. N. Pfeiffer, K. W. West, and K. W. Baldwin, Microwave spectroscopy of the low-filling-factor bilayer electron solid in a wide quantum well, Nature Communications 6, 7071 (2015).
  * Zondiner  _et al._ [2020] U. Zondiner, A. Rozen, D. Rodan-Legrain, Y. Cao, R. Queiroz, T. Taniguchi, K. Watanabe, Y. Oreg, F. von Oppen, A. Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani, Cascade of phase transitions and Dirac revivals in magic-angle graphene, Nature 582, 203 (2020).
  * Rozen  _et al._ [2021] A. Rozen, J. M. Park, U. Zondiner, Y. Cao, D. Rodan-Legrain, T. Taniguchi, K. Watanabe, Y. Oreg, A. Stern, E. Berg, P. Jarillo-Herrero, and S. Ilani, Entropic evidence for a Pomeranchuk effect in magic-angle graphene, Nature 592, 214 (2021).
  * Fogler [2001] M. M. Fogler, Stripe and bubble phases in quantum hall systems, in  _High Magnetic Fields: Applications in Condensed Matter Physics and Spectroscopy_, edited by C. Berthier, L. P. Lévy, and G. Martinez (Springer Berlin Heidelberg, Berlin, Heidelberg, 2001) pp. 98–138.
  * Lorenzana  _et al._ [2001a] J. Lorenzana, C. Castellani, and C. D. Castro, Phase separation frustrated by the long-range coulomb interaction. I. theory, Phys. Rev. B 64, 235127 (2001a).
  * Lorenzana  _et al._ [2001b] J. Lorenzana, C. Castellani, and C. Di Castro, Phase separation frustrated by the long-range coulomb interaction. II. applications, Phys. Rev. B 64, 235128 (2001b).
  * Ortix  _et al._ [2008] C. Ortix, J. Lorenzana, and C. Di Castro, Coulomb-frustrated phase separation phase diagram in systems with short-range negative compressibility, Phys. Rev. Lett. 100, 246402 (2008).
  * Ortix  _et al._ [2006] C. Ortix, J. Lorenzana, and C. Di Castro, Frustrated phase separation in two-dimensional charged systems, Phys. Rev. B 73, 245117 (2006).
  * Ortix  _et al._ [2007] C. Ortix, J. Lorenzana, M. Beccaria, and C. Di Castro, Screening effects in coulomb-frustrated phase separation, Phys. Rev. B 75, 195107 (2007).
  * Ortix  _et al._ [2009] C. Ortix, J. Lorenzana, and C. Di Castro, Universality classes for coulomb frustrated phase separation, Physica B: Condensed Matter 404, 499 (2009).
  * Drummond and Needs [2009] N. D. Drummond and R. J. Needs, Phase diagram of the low-density two-dimensional homogeneous electron gas, Phys. Rev. Lett. 102, 126402 (2009).
  * Falson  _et al._ [2022] J. Falson, I. Sodemann, B. Skinner, D. Tabrea, Y. Kozuka, A. Tsukazaki, M. Kawasaki, K. von Klitzing, and J. H. Smet, Competing correlated states around the zero-field Wigner crystallization transition of electrons in two dimensions, Nature Materials 21, 311 (2022).
  * Hossain  _et al._ [2020] M. S. Hossain, M. K. Ma, K. A. V. Rosales, Y. J. Chung, L. N. Pfeiffer, K. W. West, K. W. Baldwin, and M. Shayegan, Observation of spontaneous ferromagnetism in a two-dimensional electron system, Proceedings of the National Academy of Sciences 117, 32244 (2020).
  * Yoon  _et al._ [1999] J. Yoon, C. C. Li, D. Shahar, D. C. Tsui, and M. Shayegan, Wigner crystallization and metal-insulator transition of two-dimensional holes in gaas at $\mathit{B}\phantom{\rule{0.0pt}{0.0pt}}=\phantom{\rule{0.0pt}{0.0pt}}0$, Phys. Rev. Lett. 82, 1744 (1999).
  * Knighton  _et al._ [2018] T. Knighton, Z. Wu, J. Huang, A. Serafin, J. S. Xia, L. N. Pfeiffer, and K. W. West, Evidence of two-stage melting of wigner solids, Phys. Rev. B 97, 085135 (2018).
  * Shklovskii and Efros [1972] B. I. Shklovskii and A. L. Efros, Completely compensated crystalline semiconductor as a model of an amorphous semiconductor, Sov. Phys.-JETP 35, 610 (1972).
  * Ando  _et al._ [1982] T. Ando, A. B. Fowler, and F. Stern, Electronic properties of two-dimensional systems, Rev. Mod. Phys. 54, 437 (1982).
  * Shklovskii [2007] B. I. Shklovskii, Simple model of coulomb disorder and screening in graphene, Phys. Rev. B 76, 233411 (2007).
  * Adam  _et al._ [2007] S. Adam, E. H. Hwang, V. M. Galitski, and S. Das Sarma, A self-consistent theory for graphene transport, Proceedings of the National Academy of Sciences 104, 18392 (2007).
  * Skinner and Shklovskii [2013] B. Skinner and B. I. Shklovskii, Theory of the random potential and conductivity at the surface of a topological insulator, Phys. Rev. B 87, 075454 (2013).
  * Clark  _et al._ [2009] B. K. Clark, M. Casula, and D. M. Ceperley, Hexatic and mesoscopic phases in a 2d quantum coulomb system, Phys. Rev. Lett. 103, 055701 (2009).
  * Koulakov  _et al._ [1996] A. A. Koulakov, M. M. Fogler, and B. I. Shklovskii, Charge density wave in two-dimensional electron liquid in weak magnetic field, Phys. Rev. Lett. 76, 499 (1996).
  * Fogler  _et al._ [1996] M. M. Fogler, A. A. Koulakov, and B. I. Shklovskii, Ground state of a two-dimensional electron liquid in a weak magnetic field, Phys. Rev. B 54, 1853 (1996).
  * Ciftja and Hysi [2011] O. Ciftja and I. Hysi, The electrostatic potential of a uniformly charged disk as the source of novel mathematical identities, Applied Mathematics Letters 24, 1919 (2011).
  * Goerbig  _et al._ [2004] M. O. Goerbig, P. Lederer, and C. M. Smith, Competition between quantum-liquid and electron-solid phases in intermediate landau levels, Phys. Rev. B 69, 115327 (2004).