* Shanbhag (2019) S. Shanbhag, “Mathematical foundations of an ultra coarse-grained slip link model,” J. Chem. Phys. 151, 044903 (2019).
  * Pasquali and Scriven (2004) M. Pasquali and L. E. Scriven, “Theoretical modeling of microstructured liquids: a simple thermodynamic approach,” J. Non-Newton. Fluid Mech. 120, 101–135 (2004).
  * Stephanou, Tsimouri, and Mavrantzas (2016) P. S. Stephanou, I. C. Tsimouri, and V. G. Mavrantzas, “Flow-induced orientation and stretching of entangled polymers in the framework of nonequilibrium thermodynamics,” Macromolecules 49, 3161–3173 (2016).
  * Hulsen (1990) M. A. Hulsen, “A sufficient condition for a positive definite configuration tensor in differential models,” J. Non-Newton. Fluid Mech. 38, 93–100 (1990).
  * Robertson  _et al._ (2017) B. Robertson, R. L. Thompson, T. C. B. McLeish, and I. Robinson, “Theoretical prediction and experimental measurement of isothermal extrudate swell of monodisperse and bidisperse polystyrenes,” J. Rheol. 61, 931–945 (2017).
  * Milner (2020) S. T. Milner, “Unified entanglement scaling for flexible, semiflexible, and stiff polymer melts and solutions,” Macromolecules 53, 1314–1325 (2020).
  * Kremer and Grest (1990) K. Kremer and G. S. Grest, “Dynamics of entangled linear polymer melts: A molecular-dynamics simulation,” J. Chem. Phys. 92, 5057–5086 (1990).
  * Maklad and Poole (2021) O. Maklad and R. J. Poole, “A review of the second normal-stress difference; its importance in various flows, measurement techniques, results for various complex fluids and theoretical predictions,” J. Non-Newton. Fluid Mech. , 104522 (2021).
  * O’Connor, Alvarez, and Robbins (2018) T. C. O’Connor, N. J. Alvarez, and M. O. Robbins, “Relating chain conformations to extensional stress in entangled polymer melts,” Phys. Rev. Lett. 121, 047801 (2018).
  * Fetters, Lohse, and Colby (2007) L. J. Fetters, D. J. Lohse, and R. H. Colby, “Chain Dimensions and Entanglement Spacings,” in  _Physical Properties of Polymers Handbook_ (Springer, 2007) pp. 447–454.
  * Galvani Cunha, Olmsted, and Robbins (2022) M. A. Galvani Cunha, P. D. Olmsted, and M. O. Robbins, “Probing the nonequilibrium dynamics of stress, orientation, and entanglements in polymer melts with orthogonal interrupted shear simulations,” J. Rheol. 66, 619–637 (2022).
  * Sanchez-Reyes and Archer (2002) J. Sanchez-Reyes and L. A. Archer, “Step shear dynamics of entangled polymer liquids,” Macromolecules 35, 5194–5202 (2002).
  * Stratton and Butcher (1973) R. A. Stratton and A. F. Butcher, “Stress relaxation upon cessation of steady flow and the overshoot effect of polymer solutions,” Journal of Polymer Science: Polymer Physics Edition 11, 1747–1758 (1973).
  * Roy and Roland (2013) D. Roy and C. M. Roland, “Reentanglement kinetics in polyisobutylene,” Macromolecules 46, 9403–9408 (2013).
  * Ianniruberto and Marrucci (2014c) G. Ianniruberto and G. Marrucci, “Do repeated shear startup runs of polymeric liquids reveal structural changes?” ACS Macro Lett. 3, 552–555 (2014c).
  * Brini  _et al._ (2013) E. Brini, E. A. Algaer, P. Ganguly, C. Li, F. Rodríguez-Ropero, and N. F. A. van der Vegt, “Systematic coarse-graining methods for soft matter simulations–a review,” Soft Matter 9, 2108–2119 (2013).
  * Johnson and Phelan Jr (2021) L. C. Johnson and F. R. Phelan Jr, “Dynamically consistent coarse-grain simulation model of chemically specific polymer melts via friction parameterization,” J. Chem. Phys. 154, 084114 (2021).
  * Oldroyd (1950) J. G. Oldroyd, “On the formulation of rheological equations of state,” Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences 200, 523–541 (1950).
  * Larson (1988) R. G. Larson, “Constitutive equations for polymer melts and solutions,” (1988).
  * Snoeijer  _et al._ (2020) J. H. Snoeijer, A. Pandey, M. A. Herrada, and J. Eggers, “The relationship between viscoelasticity and elasticity,” Proc. R. Soc. A 476, 20200419 (2020).
  * Milner (1993) S. T. Milner, “Dynamical theory of concentration fluctuations in polymer solutions under shear,” Phy. Rev. E 48, 3674 (1993).


# Design of a new SPS injection system via numerical optimisation

E. Waagaard also at Uppsala University, Sweden M. J. Barnes  W. Bartmann  J.
Borburgh  L. S. Ducimetière  T. Kramer  T. Stadlbauer
P. Trubacova  F. M. Velotti  CERN  Geneva  Switzerland

###### Abstract

The Super Proton Synchrotron (SPS) injection system plays a fundamental role
to preserve the quality of injected high-brightness beams for the Large Hadron
Collider (LHC) physics program and to maintain the maximum storable intensity.
The present system is the result of years of upgrades and patches of a system
not conceived for such intensities and beam qualities. In this study, we
propose the design of a completely new injection system for the SPS using
multi-level numerical optimisation, including realistic hardware assumptions.
We also present how this hierarchical optimisation framework can be adapted to
other situations for optimal accelerator system design.

##  1 Introduction

The SPS is the last injector of the LHC, hence the second-largest accelerator
of the CERN accelerator complex [1]. One of the main duties of the SPS is to
provide high-brightness beams of high quality to the LHC physics program, as
well as high-intensity beams to the fixed-target experiments of the North
experimental Area (NA). The SPS injection system plays a fundamental role to
the SPS beam quality, not only to preserve the quality of the injected beam,
but also to maintain the maximum storable intensity in the ring. In
particular, the upcoming High-Luminosity LHC [2] requirements where the
present injection system will act as a bottleneck for higher intensities
further stresses the need for a new robust injection solution. Presently, the
beam injected into the SPS is delivered from the Proton Synchrotron (PS) via
the TT10 transfer line. The injection occurs in the horizontal plane in a
single turn. The schematic view of the main elements is illustrated in Fig. 1.

Figure 1:  Illustration of present TT10/SPS injection.

There are four septa (denoted MSI) - insertion magnets that define the
borderline between the ring and the transfer line (TL) - in the present
injection system. Four kicker magnets - fast-pulsed dipole magnets - give the
final deflection to the injected beam. There are two main types of kicker
magnets: MKP-S and MKP-L. A cross-section of the MKP-L kicker vacuum tank is
shown in Fig. 2. Due to difference in aperture, the MKP-L is more susceptible
to beam-induced heating via beam coupling impedance. Thus, the new injection
kickers should be designed also accounting for these constraints. In case of
injection failure, a dedicated beam dump denoted TBSJ is used.

Figure 2:  Cross-section of MKP-L kicker tank.

##  2 Injection simulations

For the studies and simulations, a model of the SPS injection system was
developed using MAD-X [3], presented in detail in [4]. To simplify the
integration with other frameworks for optimisation, we used the python wrapper
cpymad [5]. The three main optics used during normal SPS operations are Q20,
Q26 and SFTPRO [6]. Q20 is used for beam delivery to LHC, whereas Q26 is the
previous optics version for LHC but is now used together with Q20 for high-
intensity bunched beams of LHC type. On the other hand, SFTPRO is used for
slow extraction for the fixed-target experiments of NA.

##  3 Optimisation parameters

It is crucial to minimise errors at injection to avoid injection oscillations
in the SPS, which ultimately can lead to filamentation and emittance blow-up,
reducing luminosity and beam quality delivered to the experiments. Thus, we
strive to minimise horizontal displacement at injection with zero-valued
position $x$ and transverse momenta $x^{\prime}$ for a stable beam. A safe
beam dump is also required - when the MKPs are powered off, the beam has to
hit the TBSJ. We also need to avoid any beam losses, which occur if the beam
touches the mechanical aperture. For this purpose, the minimum acceptance
$A_{x\textrm{, min}}$ (which is the number of standard deviations of the beam
size that fits into the available space between the orbit and the mechanical
aperture) is used as a figure of merit:

$A_{x\textrm{,
min}}=\textrm{min}\bigg{(}\frac{\textrm{aper}_{x}(s)-|x(s)|}{\sigma_{x}(s)}\bigg{)},$
(1)

where $\textrm{aper}_{x}(s)$ is the horizontal mechanical aperture of the
machine, with $y$ replacing $x$ for the vertical case. For the horizontal
plane, the beam size ${\sigma_{x}(s)}$ is defined as

$\sigma_{x}(s)=\sqrt{\beta_{x}(s)\varepsilon_{x}+(D_{x}(s)\Delta p/p)^{2}},$
(2)

where $D_{x}$ is the dispersion function, $\varepsilon_{x}$ is the geometric
emittance, $\beta_{x}$ is the beta function, and $\Delta p/p$ is the momentum
spread. MAD-X simulations of the horizontal and vertical beam envelopes for
the present SPS injection system with corresponding $A_{x\textrm{, min}}$
values are shown in Fig. 3.

Figure 3:  MAD-X simulations of the current horizontal (left) and vertical
(right) envelope injection for Q20 and SFTPRO optics, with the minimum
acceptance $A_{\textrm{min}}$ for each case and at which longitudinal position
$s$ it occurs.

The strict requirements on a new injection system listed above - a stable
preserved beam that can be safely dumped, applicable to all the optics - are
primarily fulfilled by varying magnet strengths and positions of the septa and
the kickers. At the lower level, these parameters constitute the degrees of
freedom (DOF), or actors. Once these low-level requirements are fulfilled, we
also strive to deploy minimal magnet resources: maximum magnetic field $B$
required and integrated field $\int Bdl$, for facilitated magnet design. For
this reason, we can also at the higher level vary the number $n_{1}$ of MKP
magnet modules per tank and the number $n_{2}$ of MKP vacuum tanks,
illustrated in Fig. 4, to find the ideal kicker configuration for minimal
magnet resources. On the other hand, parts of the septa from the PS Booster
and parts from the present SPS injection system can be reused to form a new
movable septum conductor (also called blade), whose second half is 20 mm thick
instead of today’s 40 mm for the whole blade. The reduced thickness of the
proposed MSI blades allows for more space for the beam envelope. Thus, this
new MSI blade remains a fixed parameter.

Figure 4:  Illustration of the MKP magnet tank.

##  4 Hierarchical optimisation

The injection system requirements and actors spanning over multiple levels are
incorporated in a hierarchical optimisation framework, where the low-level
optimisation ensures a correct trajectory at injection, and the outer high-
level optimisation minimises the magnet resources. An illustration of the flow
across different hierarchy levels is shown in Fig. 5. Numerical optimisation
and machine learning are used in accelerator physics, and some attempts have
also been made to exploit genetic algorithms (GA) for multi-objective
optimisation to construct transfer lines [7], alas without a proposed final
solution.

Figure 5:  Illustration of the hierarchical optimiser.

The meta-optimiser (a wrapper function) first creates a global environment
object with a given kicker magnet configuration from the tuple ($n_{1}$,
$n_{2}$). The environment contains all the attributes and a step method, which
takes the tuple as input parameter and outputs the high-level objective
function to minimise, defined as

$\textrm{objective function}=\frac{1}{a_{1}}\bigg{(}\int
Bdl\bigg{)}^{2}+\frac{1}{a_{2}}B^{2},$ (3)

where $a_{i}$ are normalisation constants. The global environment creates low-
level environments for all the three optics (Q20, Q26 and SFTPRO). Each of
these low-level environments contain the respective simulation process in
MADX, and also a step method that takes magnet strengths and positions of the
MSI and MKP, outputting the low-level objective function

$\displaystyle\begin{split}&\textrm{objective
function}=\frac{1}{a_{3}}x^{2}+\frac{1}{a_{4}}x^{\prime 2},\\\
&\textrm{penalty}=\textrm{ max}(0,A_{x\textrm{, min}}-A_{x})^{2}\end{split}$
(4)

to which a penalty is added: the square of the excess of the present
acceptance $A_{x}$, with $A_{x\textrm{, min}}$ defined in Eq. (1). Thus, beam
losses or an improperly dumped beam trajectory are heavily penalised. The low-
level objective function with its penalties remains general, which can be
replaced by or coupled with the vertical plane if desired.

##  5 Optimisation methods

The low-level objective function is optimised for each optics with the
adaptive Nelder-Mead algorithm, implemented in scipy [8]. The global
environment then selects the magnet positions from the optics that requires
the highest $\int\leavevmode\nobreak\ Bdl$, fixes these magnet positions and
re-optimises their strengths - this approach guarantees that we first satisfy
the case requiring the highest integrated field. In our study, the magnet
strengths for the other optics could always be re-optimised with these new
magnet positions.

At the higher level, the optimiser evaluates the high-level objective function
of the global environment and finds the ideal pair ($n_{1}$, $n_{2}$), using a
GA with the library pymoo [9]. GAs are inspired by the process of natural
selection, where a starting population of initial settings is sampled and the
(high-level) objective functions are evaluated. The fittest candidates survive
and are selected for the crossover, combining into offspring where some
mutations may occur according to pre-defined probabilities. Even with a
starting population size of 10 randomly sampled pairs of settings, the
algorithm converges already after 10-15 iterations to similar candidate
solutions.

##  6 Results

The results produced by the hierarchical optimisation framework are presented
in Fig. 6, using the new thinner septum blades. The combination ($n_{1}=3$,
$n_{2}=5$) yields the lowest value of the high- and low-level objective
function values, and higher $A_{\textrm{x. min}}$ than today.

Figure 6:  Optimisation results for Q20 and SFTPRO optics (left - injected
beam, right - dumped beam).

Although $\int Bdl$ for Q20 and Q26 would increase by 5-6%, the injected un-
kicked beam reaches the TBSJ without the need of an extra dipole as used today
- a great simplification for future operation.

A close-up of the MKP tanks in the present and proposed injection system is
shown in Fig. 7. The new configuration only contains MKP-S type kickers,
leading to less beam-induced heating and more uniform design and spare policy.
The increased intermediate spacing for valves allows for independent vacuum
sectors. The last two MKP tanks have been shifted by 1.9 cm horizontally in
the new configuration to allow for more space to the dumped beam. The main
drawback of this configuration is the fact that the external side of the last
MKP tank is exposed to direct beam impact in case of self-triggering of the
MKP while beam is circulating. More studies are needed to evaluate alternative
solutions and to assess the risk associated with direct impact on the MKP,
possibly with local shielding.

Figure 7:  Present and proposed MKP configuration.

##  7 Conclusions

The current SPS injection system has seen a radical evolution since its first
design, and the characteristics of the injected beam have changed
significantly. A full refactoring may be needed to solve long-standing issues,
like the different kicker types and the bulky MSI. We presented a hierarchical
optimisation framework that gradually constructs a new injection sequence and
generates ideal candidate solutions for all SPS optics - finding the ideal
balance between minimal magnet resources and dumped beam trajectory, while
considering realistic hardware. The developed framework is rather generic and
provides a basis for designing other accelerator systems, even with higher
complexity. We stress that although the SPS injection is a specific case, the
high-level discrete design of magnet modules and configurations can in
principle be extended to almost any accelerator sequence, as described in a
pseudo-algorithm in [4].

## References

  * [1] P. Collier et al., SPS as injector for LHC: Conceptual design report (1997). Tech. Rep. CERN-SL-97-07 DI
  * [2] G. Apollinari et al., High luminosity LHC project description (2014). No. CERN-ACC-2014-0321.
  * [3] MADX Official Website. Available at: http://madx.web.cern.ch/madx/
  * [4] E. Waagaard. Designing a New SPS Injection System With Numerical Optimisation (2022). Available at http://urn.kb.se/resolve?urn=urn:nbn:se:uu:diva-468033.
  * [5] CPyMad GitHub repository, available at https://hibtc.github.io/cpymad/
  * [6] SPS Optics Repository, available at https://acc-models.web.cern.ch/acc-models/sps/.
  * [7] Y. Dutheil, F. Velotti, and B. Goddard. The blind linemaker: Genetic algorithms for multi-objective optimisation (2020). Available at https://indico.cern.ch/event/964060/.
  * [8] SciPy Official Website, available at https://scipy.org/.
  * [9] Pymoo Official Website, available at https://pymoo.org/.


# The effect of combined roundness and polydispersity on the phase behavior of
hard-rectangle fluids

Yuri Martínez-Ratón  Grupo Interdisciplinar de Sistemas Complejos (GISC),
Departamento de Matemáticas, Escuela Politécnica Superior, Universidad Carlos
III de Madrid, Avenida de la Universidad 30, E-28911, Leganés, Madrid, Spain
<EMAIL_ADDRESS>Enrique Velasco  Departamento de Física Teórica de la
Materia Condensada, Instituto de Física de la Materia Condensada (IFIMAC) and
Instituto de Ciencia de Materiales Nicolás Cabrera, Universidad Autónoma de
Madrid, E-28049, Madrid, Spain<EMAIL_ADDRESS>
###### Abstract

We introduce a model for a fluid of polydisperse rounded hard rectangles where
the length and width of the rectangular core are fixed, while the roundness is
taken into account by the convex envelope of a disk displaced along the
perimeter of the core. The diameter of the disk has a continuous
polydispersity described by a Schultz distribution function. We implemented
the scaled particle theory for this model with the aim to studying: (i) the
effect of roundness on the phase behavior of the one-component hard-rectangle
fluid, and (ii) how polydispersity affects phase transitions between
isotropic, nematic and tetratic phases. We found that roundness greatly
affects the tetratic phase, whose region of stability in the phase diagram
strongly decreases as the roundness parameter is increased. Also the interval
of aspect ratios where the tetratic-nematic and isotropic-nematic phase
transitions are of first order considerably reduces with roundness, both
transitions becoming weaker. Polydispersity induces strong fractionation
between the coexisting phases, with the nematic phase enriched in particles of
lower roundness. Finally, for high enough polydispersity and certain mean
aspect ratios, the isotropic-to-nematic transition can change from second (for
the one-component fluid) to first order. We also found a packing-fraction
inversion phenomenon for large polydispersities: the coexisting isotropic
phase has a higher packing fraction than the nematic.

##  I Introduction

The study of entropic phase transitions in two dimensional liquid crystals is
nowadays an active line of research. This is not only because the study of
orientational transitions in monolayers of anisotropic colloids or molecules
(adsorbed at surfaces, interfaces or forming membranes) continues to be an
interesting research topic review1 ; review2 ; Dogic . In addition, novel
techniques to synthesize hard-core interacting micro-prisms of any cross-
sectional geometry, subsequently adsorbed to form monolayers, have been
applied to create a plethora of effectively two-dimensional fluids of Brownian
particles with several shapes Zhao5 ; Zhao ; Zhao0 ; Zhao1 ; Zhao2 ; Zhao3 ;
Zhao4 . Examples are: squares Zhao5 , rectangles Zhao , rhombuses Zhao0 ,
triangles Zhao1 , hexagons Zhao2 , and kites Zhao3 . These two-dimensional
fluids have been a source of fascinating phenomena including the appearance of
exotic symmetries and chirality in the orientational and spatial ordering of
particles. Other exotic shapes such as circular arcs Cinacchi1 ; Cinacchi2
have also been studied via Monte Carlo (MC) simulations, showing interesting
self-assembling properties directly related to the presence of "entropic
bonding", a concept recently introduced in Ref. Glotzer . The presence of
liquid-crystal and crystalline ordering was also investigated in MC
simulations of two-dimensional hard regular polygons, where the number of
edge-lengths play a crucial role in the symmetries of the stable liquid
crystal and crystal phases Glotzer2 .

A paradigmatic example of a two-dimensional liquid-crystal which exhibits
isotropic (I), nematic (N) and tetratic (T) orientational ordering is the
fluid of hard rectangles (HR) and its hard square limit. This fluid was
extensively studied using the density functional theory (DFT) Schlaken ; MR1 ;
three-body , via MC simulations Frenkel ; Donev ; Torres , and even by
experiments consisting of vertically shaken monolayers of granular particles
of rectangular or square cross-sections Narayan ; Dani ; Miguel ; Menon .
While particles in the I phase are not orientationally ordered, in the N or T
phases the main particle axes (parallel to the major edge-length) orient along
one or two equivalent directors, respectively. The orientational distribution
function $h(\phi)$, i.e. the probability density of particles axes to align at
an angle $\phi$ with respect to the director, has twofold,
$h(\phi)=h(\phi+\pi)$ (N phase), or fourfold, $h(\phi)=h(\phi+\pi/2)$ (T
phase) symmetries. The T phase in the HR fluid has been shown to be stable
only for aspect ratios varying from 1 (the hard-square limit, with the T phase
as the only possible liquid-crystal phase) to a particular value $\kappa^{*}$.
The value of $\kappa^{*}$ predicted by scaled particle theory (SPT) MR1 , a
version of DFT, and a more sophisticated DFT based on the second and third
virial coefficients three-body , are 2.21 and 3.23, respectively. Recent
simulations have shown that $\kappa^{*}\approx 5$ Torres  while experiments on
monolayers of quasi-two-dimensional granular cylinders indicate the presence
of T correlations for aspect ratios as large as $\kappa_{c}\approx 7$ Dani ;
Miguel . However, some care should be taken to compare the results from
experiments on non-equilibrium dissipative granular rods to those obtained by
theoretical models based on equilibrium statistical mechanics. As we have
shown recently Miguel , energy dissipation in vertically-shaken granular
monolayers strongly promotes particle clusterization, resulting in a high
proportion of square-like clusters made of rectangles joined side by side,
which in turn induces the formation of T textures as stationary states.

On the other hand, confined monolayers of cylinders respond to geometrical
frustration much in the same way as equilibrium liquid crystals, i.e. by
creating topological defects that restore the global symmetry of the system.
The number and topological charge of these defects seem to follow the rules of
topology. Also, when particles are confined in annular geometries, a
complicated pattern arises in the orientational-ordering field, with domain
walls that separate regions of smectic and T ordering and additional
topological defects Ariel . While topology predicts that no defects should be
excited in this case, the small size of the system compared with the particle
length probably explains the formation of this complex pattern, although non-
equilibrium effects cannot be discarded as an origin. Colloidal monolayers
also exhibit the presence of T-like disclination defects in the smectic
textures when confined inside cavities of different shapes Lowen1 ; Lowen2 .
These similarities between dissipative and equilibrium monolayers point to the
preponderant role of entropic interactions as the main mechanism dictating the
symmetries of both systems when frustrated by confinement.

The T phase can be stabilized by other geometrical shapes such as rhombuses
Torres  and kites Zhao3 ; MR4  of particular shapes and ratios between their
characteristic lengths. Indeed its stability region in the phase diagram seems
to be very sensitive to these ratios and, what can be more important, to the
roundness of the particle corners. It was recently shown by MC simulations
that a fluid of hard rounded squares does not exhibit a T phase for high
enough roundness of the corners, with the I phase directly undergoing a
transition to a crystalline phase Escobedo . This result explained why the T
phase was not found in recent experiments on rounded squares Zhao5 : it is
certainly difficult to design an experimental procedure to obtain
microparticles with perfect right corners. There are two ways to implement the
presence of nonzero curvature in the particle boundaries: (i) to take into
account the change of curvature by defining the particle as a superellipse
Torres ; Varga1 ; Varga2  with an exponential parameter ranging from 1 (the
rhomboidal shape), 2 (the elliptical shape), and finally the infinite limit
(for the case of rectangles), or (ii) to consider a fixed core defined by
straight lines, adding the convex envelope that results when a disk of some
particular diameter slides along the boundary of this core Escobedo .

Here we will use the second recipe, with a fixed core of rectangular shape,
defining in such a way a fluid of hard rounded rectangles (HRR). Two main
studies have been carried out. In the first we study the effect of roundness
(measured through a roundness parameter) on the stability of the T phase as
compared with the HR fluid. Using scaled particle theory (SPT) Cotter ; Lasher
; Barboy , several phase diagrams for different values of the roundness
parameters have been calculated, which allowed to trace out the stability
boundaries of the different phases and their changes as a function of the
roundness parameter. The second study deals with the effect of polydispersity
in the phase behavior of HRR. This point is motivated by the fact that some
polydispersity in sizes/shapes is always present in the experimental systems.
As will be shown later, the main effects of particle roundness on the one-
component fluid of HRR are: (i) the strong destabilization of the T phase; the
stability region of the T phase in the packing fraction-aspect ratio plane is
considerably shrinked as roundness increases, (ii) the interval of aspect
ratios where the I-N and T-N transitions are of first order is strongly
reduced, and as a consequence both transitions become weaker. As regards the
effect of polydispersity we find that, for certain aspect ratios and high
enough mean roundness and polydispersity coefficient, the I-N transition for
certain mean aspect ratios becomes of first order despite being of second
order in the one-component fluid. Also when the fluid exhibits a first-order
I-N or T-N transition, and for high enough polydispersity, the coexisting
phases exhibit a packing fraction inversion due to the fractionation effect:
the coexisting I or T phases are enriched in particles of higher roundness (or
lower aspect ratios), while the N phase is more populated by species of higher
aspect ratios. As a consequence, the coexisting I or T phases can have a lower
packing fraction as compared to that of the coexisting N. We should mention
that previous MC simulation studies on monodisperse hard rods in 2D showed
that a quasilong-range ordered N phase exhibits a transition to the I phase
via a Kosterlitz-Thouless disclination unbinding type mechanism rather than
being of first order Frenkel2 ; Dijkstra . However recent studies have shown
that, for particular types of particle interactions, the I-N transition
becomes of first order in 2D Enter ; Vink . Finally, recent experiments found
that quasi-monolayers of magnetic nanorods confined between adjacent layers of
a lamellar phase exhibit a first order I-N transition Constantin .

The theoretical DFT study of continuous polydisperse fluids of anisotropic
particles represents a challenge because the density profile depends not only
on the particle orientational degrees of freedom. It also incorporates the
distribution of the polydisperse variable, which complicates the numerical
procedure necessary to calculate phase coexistence. To deal with this problem,
some simplifications were made in the past to study the effect of
polydipersity on the phase behavior of freely-rotating hard polydisperse rods.
One of these simplifications involves using the Onsager-DFT of hard
spherocylinders in the hard-needle limit, and implementing the spherical
harmonics expansion of the excluded volume up to second order, together with
the use of the moment theory Sollich1a ; Sollich1b  to render the calculations
feasible Sollich2 . An alternative approach is to discretize the orientational
degrees of freedom, as in the Zwanzig approximation, and use the Fundamental
Measure DFT for hard boardlike particles which correctly describes not only
two-body, but also three-body correlations Cuesta . As we show here the
present model has the advantage that the orientational degrees of freedom and
the polydisperse variable (the diameter of the disk causing particle
roundness) are decoupled, making the theoretical treatment of the polydisperse
fluid as easy as its one-component counterpart. This property can be used in
future developments to study the combined effect of confinement and
polydispersity on the structural properties of a two-dimensional liquid-
crystal fluid. The article is organized as follows. Sec. II is devoted to the
definition of the model and the presentation of the theoretical tools used for
the calculation of phase diagrams. Special attention is paid to the
coexistence calculation formalism (Sec. II.1), the definition of the
polydisperse distribution function used in the study (Sec. II.2), and the
implementation of the bifurcation analysis to calculate the second-order phase
transition curves (Sec. II.3). The results are divided in two parts: In Sec.
III we present the results for the one-component fluid (zeroth
polydispersity), while in Sec. IV we describe the results regarding the effect
of polydispersity on the phase behavior of HRR. In Sec. V we describe an
approximate procedure to account for the effect of roundness on the
instability of the T phase with respect to crystallization. Finally some
conclusions and discussions are summarized in Sec. VI. We relegate to Sec. A
the details for the numerical calculations of shadow and cloud coexistence
curves in the polydisperse HRR fluid.

##  II Model and Theory

Figure 1: Sketches of rounded rectangles where the fixed core lengths $L$ and
$\sigma$ and the polydisperse disk length $l$ of particle are indicated. Note
that the limiting case $l=0$ corresponds to a rectangular particle, while for
$l\neq 0$ and $\sigma=0$ one obtains a discorectangle.

Our model consists of a polydisperse mixture of hard rounded rectangles (HRR).
A particle is defined by a fixed rectangular core of length $L$ and width
$\sigma$ ($L\geq\sigma$) plus a polydisperse coating obtained by making the
center of a disk of diameter $l$ to slide over the perimeter of the
rectangular core (see Fig. 1 for a sketch of the particle geometry). While the
core sizes $L$ and $\sigma$ are fixed, the diameter $l$ is a polydisperse
variable with a value ranging from 0 to $\infty$. Note that the cases $l=0$
and $\sigma=0$ constitute the limiting cases of hard rectangular and
discorectangular particles, respectively. In the following we define the main
axis of a particle to be parallel to its length $L$.

Let us consider two such particles with different values ($l$ and
$l^{\prime}$) of the disk diameters. The excluded area between these
particles, as a function of their relative angle
$\Delta\phi=\phi-\phi^{\prime}$, can be computed as

$\displaystyle A_{\rm
excl}(\Delta\phi,l,l^{\prime})=\left(L^{2}+\sigma^{2}\right)|\sin\Delta\phi|+2L\sigma|\cos\Delta\phi|$
$\displaystyle+\left(L+\sigma\right)\left(l+l^{\prime}\right)+\frac{\pi}{2}ll^{\prime}+a(l)+a(l^{\prime}),$
(1)

where the particle area, $a(l)$, is

$\displaystyle a(l)=L\sigma+(L+\sigma)l+\frac{\pi}{4}l^{2}.$ (2)

From the excluded area we define the geometric function from which the scaled
particle theory is constructed:

$\displaystyle A_{\rm spt}(\Delta\phi,l,l^{\prime})=\frac{1}{2}\left[A_{\rm
excl}(\Delta\phi,l,l^{\prime})-a(l)-a(l^{\prime})\right].$ (3)

The polydisperse mixture of HRRs is characterized by its number density
distribution function, $\rho(l,\phi)$, a probability density in the variables
$l$ and $\phi$. This function describes the probability to find a particle
with disk size $l$ oriented with an angle $\phi$ with respect to a fixed
reference frame. Thus we have
$\displaystyle{\int_{0}^{\infty}dl\int_{0}^{2\pi}d\phi\rho(l,\phi)=\frac{N}{A}}$,
with $N$ and $A$ the total number of particles and the area of the system.
From $\rho(l,\phi)$ we can define its $i$th moment $m_{i}(\phi)$ with respect
to $l$ and the integral of this moment with respect to $\phi$ as

$\displaystyle m_{i}(\phi)\equiv\int_{0}^{\infty}dl\ l^{i}\rho(l,\phi),\quad
m_{i}\equiv\int_{0}^{2\pi}d\phi\ m_{i}(\phi).$ (4)

Note that $m_{0}=N/A$, the total number density. To facilitate numerical
calculations we will use the Fourier expansion of the function $m_{0}(\phi)$,

$\displaystyle m_{0}(\phi)=\frac{1}{2\pi}\left[m_{0}+\sum_{k\geq
1}m_{0}^{(k)}\cos(2k\phi)\right],$ (5)

with $\\{m_{0}^{(k)}\\}$ the Fourier amplitudes. From $m_{0}(\phi)$ the
orientational distribution function is simply $h(\phi)=m_{0}(\phi)/m_{0}$,
while the order parameters describing the orientational ordering are

$\displaystyle Q_{2n}\equiv\int_{0}^{2\pi}d\phi
h(\phi)\cos(2n\phi)=\frac{m_{0}^{(n)}}{2m_{0}},\quad n=1,2$ (6)

For uniaxial N orientational symmetry we have $Q_{2}\neq 0$, while T symmetry
is characterized by $Q_{2}=0$ and $Q_{4}\neq 0$.

The double average of the function $A_{\rm
spt}(\phi-\phi^{\prime},l,l^{\prime})$ with respect to $\rho(l,\phi)$ and
$\rho(l^{\prime},\phi^{\prime})$ gives

$\displaystyle\langle\langle A_{\rm
spt}(\Delta\phi,l,l^{\prime})\rangle\rangle_{\rho(l,\phi)}\equiv\int_{0}^{\infty}dl\int_{0}^{\infty}dl^{\prime}\int_{0}^{2\pi}d\phi$
$\displaystyle\times\int_{0}^{2\pi}d\phi^{\prime}\rho(l,\phi)\rho(l,\phi^{\prime})A_{\rm
spt}(\phi-\phi^{\prime},l,l^{\prime})$
$\displaystyle=\frac{\left(L^{2}+\sigma^{2}\right)}{2}\langle\langle|\sin(\Delta\phi)|\rangle\rangle_{m_{0}(\phi)}$
$\displaystyle+L\sigma\langle\langle|\cos(\Delta\phi)|\rangle\rangle_{m_{0}(\phi)}+(L+\sigma)m_{0}m_{1}+\frac{\pi}{4}m_{1}^{2},$
(7)

where we have used the shorthand notation

$\displaystyle\langle\langle
g(\Delta\phi)\rangle\rangle_{m_{0}(\phi)}=\int_{0}^{2\pi}d\phi\int_{0}^{2\pi}d\phi^{\prime}m_{0}(\phi)m_{0}(\phi^{\prime})$
$\displaystyle\times g(\phi-\phi^{\prime}).$ (8)

Inserting the Fourier expansion (5) into Eq. (7) we obtain

$\displaystyle\langle\langle A_{\rm
spt}(\Delta\phi,l,l^{\prime})\rangle\rangle_{\rho(l,\phi)}=\frac{g_{0}m_{0}^{2}-\frac{1}{2}\sum_{k\geq
1}g_{k}\left(m_{0}^{(k)}\right)^{2}}{\pi}$
$\displaystyle+(L+\sigma)m_{0}m_{1}+\frac{\pi}{4}m_{1}^{2},$ (9)

with

$\displaystyle g_{k}\equiv\frac{\left(L+(-1)^{k}\sigma\right)^{2}}{4k^{2}-1}$
(10)

Another important quantity of the polydisperse mixture is the total packing
fraction:

$\displaystyle\eta=\int_{0}^{\infty}dl\int_{0}^{2\pi}d\phi\rho(l,\phi)a(l)=m_{0}L\sigma+m_{1}\left(L+\sigma\right)$
$\displaystyle+\frac{\pi}{4}m_{2},$ (11)

In this expression we used (2) for the particle area and (4) for the
integrated moments of $\rho(l,\sigma)$. Note that $\eta$ depends not only on
the the zeroth and first moments $m_{0}$ and $m_{1}$, as does the double
average of the scaled-particle area (9), but also on the second moment
$m_{2}$.

With these definitions, the excess part of the free-energy density, according
to the SPT three-body , can be calculated (in thermal units), as

$\displaystyle\Phi_{\rm ex}[\rho(l,\phi)]\equiv\frac{\beta{\cal F}_{\rm
exc}[\rho(l,\phi)]}{A}=-m_{0}\log(1-\eta)$ $\displaystyle+\frac{\langle\langle
A_{\rm spt}(\Delta\phi,l,l^{\prime})\rangle\rangle_{\rho(l,\phi)}}{1-\eta},$
(12)

while the ideal part is, as usual,

$\displaystyle\Phi_{\rm id}[\rho(l,\phi)]\equiv\frac{\beta{\cal F}_{\rm
id}[\rho(l,\phi)]}{A}$
$\displaystyle=\int_{0}^{\infty}dl\int_{0}^{2\pi}\rho(l,\phi)\left[\log\rho(l,\phi)-1\right].$
(13)

In the above expressions, $\beta{\cal F}_{\rm id,exc}[\rho(l,\phi)]$ are the
ideal and the excess parts of the free-energy density functional, scaled with
the factor $\beta=\left(k_{B}T\right)^{-1}$ (note that thermal area inside the
logarithm of the ideal part has been dropped).

###  II.1 Coexistence calculations

Now we calculate the two-phase coexistence between a phase that occupies a
fraction $1-\epsilon$ of the total area (the cloud phase), and another phase
that occupies a vanishingly small fraction of the area, $\epsilon\ll 1$ (the
shadow phase), with coexisting density distributions $\rho_{\rm c}(l,\phi)$
and $\rho_{\rm s}(l,\phi)$, respectively.

Let us obtain the equations that govern this coexistence. Mass conservation,
expressed by the lever rule, states that the sum of the two density
distributions, integrated over the angle $\phi$ and each multiplied by its
respective area occupancy, $1-\epsilon$ or $\epsilon$, is a conserved
quantity. This is equal to the distribution function of the parent phase,
$\rho_{0}(l)\equiv\rho_{0}f(l)$, where $\rho_{0}$ is the total number density
of the system, $\rho_{0}=N/A$, while $f(l)$ is a fixed probability disk-
diameter distribution function. The lever rule is then

$\displaystyle\rho_{0}f(l)=(1-\epsilon)\int_{0}^{2\pi}d\phi\rho_{c}(l,\phi)+\epsilon\int_{0}^{2\pi}d\phi\rho_{s}(l,\phi).$
(14)

Minimizing the total free-energy density $\Phi[\rho(l,\phi)]=\Phi_{\rm
id}[\rho(l,\phi)]+\Phi_{\rm ex}[\rho(l,\phi)]$ with respect to
$\rho_{c,s}(l,\phi)$, and using the lever rule (14) and the integral
expression (7) for the averaged scaled particle area, we obtain

$\displaystyle\rho_{c,s}(l,\phi)=e^{\beta\mu_{0}(l)}e^{-c_{1}^{(c,s)}(l,\phi)},$
(15)

where the Lagrange multiplier $\beta\mu_{0}(l)$, necessary to satisfy the
constraint (14), is just the scaled chemical potential of the species with
disk-diameter value $l$. In the above we have used the notation
$c_{1}^{(\alpha)}(l,\phi)$ for the first functional derivative of the excess
free energy:

$\displaystyle c_{1}^{(\alpha)}(l,\phi)=\frac{\delta\Phi_{\rm
ex}[\rho_{\alpha}(l,\phi)]}{\delta\rho_{\alpha}(l,\phi)}$
$\displaystyle=-\log\left(1-\eta_{\alpha}\right)+\frac{1}{1-\eta_{\alpha}}\times\frac{\delta\langle\langle
A_{\rm
spt}(\Delta\phi,l,l^{\prime})\rangle\rangle_{\rho_{\alpha}(l,\phi)}}{\delta\rho_{\alpha}(l,\phi)}$
$\displaystyle+\beta p_{\alpha}a(l),$ (16)

where

$\displaystyle\beta
p_{\alpha}=\frac{m_{0}^{(\alpha)}}{1-\eta_{\alpha}}+\frac{\langle\langle
A_{\rm
spt}(\Delta\phi,l,l^{\prime})\rangle\rangle_{\rho_{\alpha}(l,\phi)}}{(1-\eta_{\alpha})^{2}},\
\alpha=c,s,$ (17)

is the pressure of the coexisting $\alpha$-phase. Using the definition (7) and
the Fourier expansion (5), we explicitly find the first functional derivative
of the scaled particle area:

$\displaystyle\frac{\delta\langle\langle A_{\rm
spt}(\Delta\phi,l,l^{\prime})\rangle\rangle_{\rho_{\alpha}(l,\phi)}}{\delta\rho_{\alpha}(l,\phi)}$
$\displaystyle=\frac{2}{\pi}\left[g_{0}m_{0}^{(\alpha)}-\sum_{k\geq
1}g_{k}m_{0}^{(k,\alpha)}\cos(2k\phi)\right]$
$\displaystyle+\left(L+\sigma\right)m_{1}^{(\alpha)}+\left[\left(L+\sigma\right)m_{0}^{(\alpha)}+\frac{\pi}{2}m_{1}^{(\alpha)}\right]l.$
(18)

The Lagrange multiplier $\mu_{0}(l)$ can be found by inserting (15) into the
lever rule (14), and taking the limit $\epsilon\to 0$, which allows us to
obtain expressions for the coexisting cloud and shadow densities:

$\displaystyle\rho_{\alpha}(l,\phi)=\rho_{0}f(l)\frac{e^{-c_{1}^{(\alpha)}(l,\phi)}}{\displaystyle\int_{0}^{2\pi}d\phi^{\prime}e^{-c_{1}^{(c)}(l,\phi^{\prime})}},\
\alpha=c,s.$ (19)

Finally, multiplying Eqns. (19) by $l^{i}\cos(2k\phi)$, and integrating over
$\phi$ and $l$, we find

$\displaystyle
m_{i}^{(k,\alpha)}=\frac{2\rho_{0}}{1+\delta_{k0}}\int_{0}^{\infty}dll^{i}f(l)$
$\displaystyle\times\frac{\displaystyle\int_{0}^{2\pi}d\phi\cos(2k\phi)e^{-c_{1}^{(\alpha)}(l,\phi)}}{\displaystyle\int_{0}^{2\pi}d\phi^{\prime}e^{-c_{1}^{(c)}(l,\phi^{\prime})}},\
\alpha=c,s;\ i=0,1,2,$ (20)

where $m_{i}^{(k,\alpha)}$ is defined as the $k$th-order Fourier coefficient
of the moment (4).

For the cloud-coexisting phase and $k=0$, we obtain
$m_{i}^{(0,c)}=m_{i}^{(c)}=\rho_{0}\int_{0}^{\infty}dlf(l)l^{i}=\rho_{0}\langle
l^{i}\rangle_{f(l)}$, which coincides with the $i$th-moment of the parent
distribution function. Moreover if the cloud-phase is I, we have
$m_{i}^{(k,c)}=0$ $\forall$ $k\geq 1$. We have solved a subset of Eqns. (20),
together with the pressure equality, $p^{(c)}=p^{(s)}$, between cloud and
shadow phases, to find the set of moments $\\{m_{i}^{(k,\alpha)}\\}$ in both
coexisting phases and at the parent number density $\rho_{0}$. As will be
shown in Sec. A, we need to solve a total number of equations less than that
in (20), which is a direct consequence of the peculiar form of the spt-area
(3).

###  II.2 The polydisperse probability parent distribution function

In the present study we use a Schultz distribution to describe the
polydispersity in $l$ in the parent phase:

$\displaystyle
f(l)=\frac{(\nu+1)^{\nu+1}}{l_{0}\Gamma(\nu+1)}\left(\frac{l}{l_{0}}\right)^{\nu}e^{-(\nu+1)l/l_{0}},$
(21)

where $l_{0}=\langle l\rangle_{f(l)}$ is the mean value. $\Gamma(x)$ is the
Gamma function. The above expression fulfills the normalization condition
$\int_{0}^{\infty}dlf(l)=1$. The parameter $\nu\in[0,\infty]$ is related to
the mean square deviation by

$\displaystyle s\equiv\sqrt{\frac{\langle l^{2}\rangle_{f(l)}}{\langle
l\rangle_{f(l)}^{2}}-1}=\frac{1}{\sqrt{\nu+1}}.$ (22)

The parameter $s\in[0,1]$ is used as a measure of polydispersity. Note that
the second moment is given in terms of $s$ by $\langle
l^{2}\rangle_{f(l)}=l_{0}^{2}(1+s^{2})$.

To measure the fractionation between the coexisting phases, we will use the
mean disk size with respect to the distribution function of the shadow phase:

$\displaystyle
f^{(s)}(l)\equiv\frac{1}{m_{0}^{(s)}}\int_{0}^{2\pi}d\phi\rho^{(s)}(l,\phi),$
(23)

and scaled with the mean size $l_{0}$ of the cloud phase. The result is

$\displaystyle\frac{\langle
l\rangle_{f^{(s)}(l)}}{l_{0}}=\frac{m_{1}^{(s)}}{m_{0}^{(s)}l_{0}}.$

Obviously this magnitude is unity for the cloud phase.

###  II.3 I-N and I-T bifurcations

The packing fraction at bifurcation between the I phase and the
orientationally ordered N or T phases can be obtained from a bifurcation
analysis (see the Appendix A). This value gives the exact location of the
second-order transition, and corresponds to the spinodal instability of the I
phase in the case of a first-order transition. Here we express the result
given by Eqn. (45) for the packing fractions at bifurcation as a function of
new parameters characterizing particle geometry. We define the mean roundness
parameter $\theta$ and the mean aspect ratio $\kappa$ of the particle as

$\displaystyle\theta\equiv\frac{l_{0}}{\sigma+l_{0}},\quad\kappa\equiv\frac{L+l_{0}}{\sigma+l_{0}}.$
(24)

We can see that $\theta=0$ for the perfect rectangular particle ($l_{0}=0$)
and $\theta\to 1$ for a very large roundness, $l_{0}\gg\sigma$. In terms of
these new variables, the I-N ($k=1$) and I-T ($k=2$) bifurcation point takes
place at packing fractions given by

$\displaystyle\eta_{1}=\left\\{1+\frac{2(\kappa-1)^{2}}{3\pi\left[\kappa-\theta^{2}\left(1-\pi(1+s^{2})/4\right)\right]}\right\\}^{-1},$
(25)
$\displaystyle\eta_{2}=\left\\{1+\frac{2(\kappa+1-2\theta)^{2}}{15\pi\left[\kappa-\theta^{2}\left(1-\pi(1+s^{2})/4\right)\right]}\right\\}^{-1}.$
(26)

The crossover aspect ratio $\kappa_{c}$ when the I-N and I-T bifurcation
curves coincide, i.e. for $\eta_{1}(\kappa)=\eta_{2}(\kappa)$, is a linear
function of $\theta$:

$\displaystyle\kappa_{c}=\kappa^{*}-(\kappa^{*}-1)\theta,\quad\kappa^{*}=\frac{3+\sqrt{5}}{2},$
(27)

where $\kappa^{*}\approx 2.618$ is the cross-over aspect ratio corresponding
to hard rectangles ($\theta=0$). For $\kappa<\kappa_{c}$ ($\kappa>\kappa_{c}$)
the stable phase above the bifurcation curve is T (N). We can see that
$\kappa_{c}$ decreases with $\theta$ indicating that the roundness
destabilizes the T phase. The packing fraction value at $\kappa_{c}$ is

$\displaystyle\eta_{c}=\eta_{1}(\kappa_{c})=\eta_{2}(\kappa_{c})$
$\displaystyle=\left\\{1+\frac{2\kappa^{*}(1-\theta)^{2}}{3\pi\left[(\kappa^{*}+\theta)(1-\theta)+\pi\theta^{2}(1+s^{2})/4\right]}\right\\}^{-1}.$
(28)

Also, for polydisperse rectangles close to the hard-square shape ($\kappa=1$),
the I phase bifurcates to the T phase at packing fraction

$\displaystyle\eta_{2}(1)=\left\\{1+\frac{8(1-\theta)^{2}}{15\pi\left[1-\theta^{2}\left(1-\pi(1+s^{2})/4\right)\right]}\right\\}^{-1}.$
(29)

Figure 2: I-N ($\eta_{1}(\kappa)$ with dashed lines) and I-T
($\eta_{2}(\kappa)$ with solid lines) bifurcation curves in the one-component
fluid of rounded hard rectangles. The solid circles indicate the location of
the (I-N)–(I-T) crossover points at $\kappa_{c}$, while solid squares
correspond to the I-T bifurcation packing fractions for hard squares. The
inset show the crossover packing fraction value
$\eta_{c}=\eta_{i}(\kappa_{c})$, $i=1,2$ (dashed line) and the packing
fraction at the I-T bifurcation, $\eta_{2}(1)$, for hard squares, both as a
function of the roundness $\theta$.

We firstly analyze the case $s=0$, the one-component fluid. In Fig. 2 we plot
the functions $\eta_{1}(\kappa)$ (dashed lines) and $\eta_{2}(\kappa)$ (solid
lines) for values of $\theta$ belonging to the set $\\{0,\ 0.15,\ 0.3,\ 0.5,\
0.8,\ 1\\}$. We can see that the N phase stabilizes at lower densities as
$\kappa>\kappa_{c}$ increases, while for $\kappa<\kappa_{c}$ the I phase
bifurcates to the T phase. Note that the function $\eta_{1}(\kappa)$ for a
fixed value of $\kappa$ is a decreasing function of $\theta$ although the
variation is rather small. This in turn means that the roundness keeps
approximately the same I-N bifurcation value. From the figure we confirm that
the region of T phase stability strongly decreases with $\theta$ and
disappears altogether for $\theta=1$ (hard discorectangles): Not only the
value of $\kappa_{c}$ decreases with $\theta$ but also the packing fraction
$\eta_{2}(\kappa)$ dramatically increases with $\theta$. In the inset the
functions $\eta_{i}(\kappa_{c})$, $i=1,2$ (dashed) and $\eta_{2}(1)$ (solid)
are plotted as a function of $\theta$, both being monotonically increasing
functions, a direct consequence of the destabilizing effect of roundness on
the T phase.

Figure 3: The difference,
$\Delta\eta(\kappa)=\eta_{i}(\kappa;s)-\eta_{i}(\kappa;0)$ ($i=1$: dashed, and
$i=2$: solid), between the bifurcation packing fractions of polydisperse and
one-component fluid for a fixed value of the roundness $\theta=0.3$ and
different values of the polydisperse coefficient as they are shown.

For a fixed roundness $\theta$ the polydispersity has the effect of increasing
the packing fractions at I-N and I-T bifurcations, which can be seen in Fig. 3
where the difference,
$\Delta\eta(\kappa)\equiv\eta_{i}(\kappa;s)-\eta_{i}(\kappa;0)$ ($i=1,2$),
between the bifurcation packing fractions of polydisperse and one-component
fluids for the selected set of polydisperse coefficients $s=\\{0.3,\ 0.5,\
0.7,\ 1\\}$ are plotted. From the figure we conclude that this effect is
rather small.

As we will show in Sec. III, the I-N transition is of first order for
$\kappa_{c}\lesssim\kappa<\kappa_{t}$, with $\kappa_{t}$ the aspect ratio
value of the I-N tricritical point, its value strongly depending on $\theta$.
For this range of $\kappa$ we should bear in mind that the packing fractions
at which the N phase begins to be stable do not coincide with the bifurcation
values calculated here. Also when polydispersity is large enough, the shadow
and cloud curves have a large coexisting gap, strongly deviating from the
bifurcation curves.

##  III The one-component fluid

Figure 4: Phase diagrams packing fraction, $\eta$, vs. aspect ratio, $\kappa$,
of the one-component HR (a) and RHR (b) fluids. The roundness parameter for
the later is $\theta=0.3$. The regions of stability of I, T and N phases are
correspondingly labeled. With solid and dashed lines we show second and first
order phase transition. For the later the coexistence gaps are shaded in grey.
The regions of stability of the I, N and T phases are labeled.

In this section we present results for the phase behavior of the one-component
fluid of HRR. The polydisperse coefficient is set as $s=0$ in the system of
Eqns. (31)-(33), which are solved numerically, together with the pressure
equality condition, to calculate the coexistence between the I or T phases and
the N phase. In the case of a second-order transition the expressions (25) and
(26) are used to compute the packing fraction at the I-N and I-T bifurcations,
respectively, or else Eqn. (47) to find the T-N bifurcation numerically. Fig.
4(a) shows the phase diagram of the HR fluid ($\theta=0$) already obtained in
Ref. MR1 , which is plotted here for the sake of comparison. In panel (b) the
phase diagram of HRR with roundness $\theta=0.3$ is shown. The main
differences between the phase diagrams are: (i) The region of stability of the
T phase of HRR shrinks considerably (the second-order I-T bifurcation moves to
higher packing fraction substantially, while the end critical point (the point
at which the I-T second-order line and the I binodal of the I-N coexistence
meet) moves to smaller aspect ratios. This result is in agreement with the
evolution of the bifurcation curves as $\theta$ increases, a point already
discussed in Sec. II.3. (ii) The range of particle aspect ratios located
between the T-N and I-N tricritical points (the left and right ends of the
continuous lines) is considerably smaller as compared to the HR fluid, i.e.
the interval in $\kappa$ for which the I-N transition is of first order
strongly reduces. Note that the aspect ratio at the I-N tricritical point is
smaller than its HR counterpart. Not only that: also the coexistence gap
(compare the grey shaded regions inside the coexisting binodals in both
panels) is much smaller, indicating a weaker first-order transition. We can
conclude that particle roundness strongly destabilizes the T phase, making the
fluid prone to exhibiting continuous phase transitions. These features become
more pronounced as the roundness parameter $\theta$ increases. In the limit
$\theta\to 1$ we obtain the phase diagram of hard discorectangles, with a
simple second-order I-N transition for any aspect ratio.

Figure 5: N ($Q_{2}$) and T ($Q_{4}$) order parameters along the I-N, I-T
coexisting binodals and also along the T-N second order transition
corresponding to a fluid of HRR with roundness $\theta=0.3$. Solid circles
indicate the positions of the T-N (left) and I-N (right) tricritical points
while the solid square is at the critical end-point where the I-T second order
line meets the T binodal of the T-N coexistence (at the left of the point) and
the I binodal of the I-N transition (at the right). See phase diagram of Fig.
4 (b).

We now proceed to describe the orientational ordering along the coexisting and
second-order curves of the phase diagram in Fig. (4)(b). In Fig. 5 the order
parameters $Q_{2}$ and $Q_{4}$ along these curves are shown. $Q_{2}$ is
different from zero between the T-N (at $\kappa=1.68$) and I-N (at
$\kappa=3.91$) tricritical points, shown with solid circles in the figure, and
exhibits a maximum around $\kappa\approx 2.17$ where the I-N transition is
strongly of first order. At the left of $\kappa\approx 1.86$ (the location of
the critical end-point where the second order I-T transition and the first-
order I-N transition meet, indicated by a solid square), the I-N transition
continues as a T-N transition with the order parameter $Q_{4}$ of the T phase
increasing as $\kappa$ decreases up, to the intersection with $Q_{4}$ of the N
phase at the T-N tricritical point. For still lower values of $\kappa$ the T-N
transition is always of second order and the T order parameter $Q_{4}$
increases along the T-N bifurcation curve up to a value of unity at
$\kappa=1$.

Figure 6: (a) Order parameters $Q_{2}$ and $Q_{4}$ of the coexisting N phase
along the I-N coexistence of HRR with fixed core length $L=2.5$ and total
width $\sigma+l=1$ (implying $\theta=l$) as a function of $\theta$. Inset:
Packing fraction difference between the I and N phases along coexistence. (b)
The same order parameters but this time along the T-N coexistence of HRR with
$L=1.7$ and same total width $\sigma+l=1$.

It is important to quantify the topological changes in the phase diagrams of
HRR when the roundness $\theta$ is changed. A possible way to achieve this is
to calculate how the aspect ratios at the tricritical and critical end-points
change as a function of $\theta$. With this information we can trace out the
boundaries where first- and second-order transitions take place. To find these
multi-critical points we implemented the following procedure: (i) we fixed the
core length $L$ and the total width of the particle to unity: $\sigma+l=1$,
implying a roundness $\theta=l$ and an aspect ratio $\kappa=L+\theta$. Next we
changed $\theta$ from an initial value where a stable I-N or T-N coexistence
exists, and move along the coexistence curves in the direction where the
transition weakens, eventually ending in the tricritical and critical end-
points. By computing the order parameters $Q_{2}$ and $Q_{4}$ as a function of
$\theta$, and extrapolating their values to zero, we approximately obtain the
locations of these points. In Fig. 6 two examples of this procedure are shown.
In panel (a) the evolution of $Q_{2}$ and $Q_{4}$ along the I-N transition is
plotted for particles with $L=2.5$, from the initial value $\theta=0$ up to
$\theta\approx 0.5$. Beyond this point the transition becomes of second order,
and the I-N tricritical point can be identified. In (b) we selected $L=1.7$
and and the initial value $\theta=0.17$. Moving to the right along the T-N
coexistence, and extrapolating $Q_{4}$ of the T phase to zero, we obtain the
value of $\theta$ corresponding to the critical end-point beyond which the T-N
transition turns into the I-N transition. Now going to the left, and
extrapolating $Q_{2}$ of the N phase to zero, we obtain the T-N tricritical
point beyond which the T-N transition becomes of second order. Repeating this
procedure for different values of the core length $L$, one can trace out the
location of three lines in the $\theta$-$\kappa$ plane: two lines correspond
to the I-N and T-N tricritical points, and the third is identified with
critical end-points. These lines are the boundaries of the regions where
first- and second-order I-N and T-N transitions can be found. The result is
shown in Fig. 7(a). As $\theta$ increases the range of aspect ratios where
first-order transitions occur considerably shrinks, disappearing in the limit
of hard discorectangles, $\theta\to 1$. In panel (b) the packing fractions
along the multi-critical curves are represented. Clearly their packing
fraction values increase dramatically with $\theta$ implying that, for high
enough values of $\theta$, the regions of first-order transitions become
unstable with respect to crystallization of the liquid-crystal fluid.

These figures also show the loss of stability of the T phase with $\theta$:
Note that in panel (a) this phase is present in the regions labeled as 2TN and
1TN, both of which considerably shrink with $\theta$. Also in panel (b) it is
apparent that the packing fraction values beyond which the T phase is stable
(the location of the critical end-points in dashed line) increases with
$\theta$, confirming the destabilizing effect of roundness on the T phase.

Figure 7: (a): Regions spanned by the roundness, $\theta$, and aspect ratio
$\kappa$ of RHR where the I-N transition is of first (1IN) or second (2IN)
order and where the T-N transition is of first (1TN) and second (2TN) order.
(b): Packing fraction, $\eta$, at I-N (solid) and T-N (dashed) tricritical
points as a function of $\theta$. With dotted line we show the value of $\eta$
corresponding to the end-critical point.

##  IV Polydisperse HRR

Figure 8: (a) Packing fractions $\eta$ vs. polydispersity coefficient $s$ for
HRR with core length and width equal to $L=3$ and $\sigma=0.6$, respectively,
and with mean disk diameter $l_{0}=1$. The resulting mean aspect ratio and
roundness are $\kappa=2.5$ and $\theta=0.625$, respectively. The shadow and
cloud coexistence curves (for both I and N) are correspondingly labeled. The
inset shows the fractionation at the I and N shadow coexisting phases measured
through the quantity $\langle l\rangle_{f^{(s)}(l)}/l_{0}$ (see text for
definition) as a function of $s$. (b) Order parameters $Q_{2}$ and $Q_{4}$ as
a function of polydispersity along the cloud and shadow coexistence curves.

In this section we study the effect of roundness polydispersity on the phase
transitions of HRR. It is first shown that, when the mean roundness $\theta$
and the polydispersity $s$ are high, the I-N transition can be of first order.
Note that this transition is of second order in the one-component fluid. First
we solved the set of Eqns. (31)-(33), together with the pressure equality of
the coexisting phases, to find the shadow and cloud I and N coexisting curves
for those values of $s$ where the I-N transition is of first order. Also, we
used Eqn. (28) which provides an analytical expression for the packing
fraction at the second-order I-N transition. We selected a particular case of
core length and width with $L=3$ and $\sigma=0.6$, while the mean disk
diameter is fixed to $l_{0}=1$. The resulting mean aspect ratio and roundness
(see Eqns. (24)) are $\kappa=2.5$ and $\theta=0.625$, respectively. The
polydisperse coefficient $s$ was varied from zero to unity, and the packing
fraction at the I-N bifurcation or two-phase coexistence was calculated.

The results are plotted in Fig. 8 (a). The transition is of second order up to
a polydispersity of $s\simeq 0.55$, which is a tricritical point. Beyond this
point the I-N transition is of first order, with the cloud-I and cloud-N
curves inside their shadow counterparts. Note also how the shadow-N (cloud-I)
curve has a lower (higher) value of $\eta$ than the cloud-I (shadow-N) curve.
This means that the I phase is always enriched in particles with larger
roundness than in the N phase, i.e. there is fractionation in the two-phase
coexistence. Since the packing fraction depends on the mean values of $l$ and
$l^{2}$ through the moments $m_{1}$ and $m_{2}$ (see Eqn. (11), the phase with
the largest roundness will have a higher packing fraction. As a result, the
orientationally ordered N phase is populated by particles with a higher mean
aspect ratio $\kappa$. This effect can be better visualized in the inset of
panel (a), where the mean disk diameter, scaled with $l_{0}$, is plotted. The
average is taken with respect to the shadow-I or shadow-N length distribution
functions $f^{(s)}(l)$. The inset confirms the enhancement of fractionation
with $s$: the mean disk diameter is an increasing (decreasing) function of $s$
along the shadow-I (shadow-N) coexisting curve. We should bear in mind that
the scaled mean diameter along the cloud coexisting curves is always equal to
unity. The orientational ordering along coexistence is shown in panel (b). The
order parameter $Q_{2}$ departs from zero at the tricritical point and follows
the usual square-root law, $Q_{2}\sim a\sqrt{s-s_{c}}$, in the neighbourhood
of $s_{c}$. By contrast, $Q_{4}$ follows a linear trend. It is interesting to
note how both, cloud-N and shadow-N coexisting phases, have similar
orientational orderings, with the cloud phase having a slightly higher values
of $Q_{2n}$ than the shadow phase, in particular for $s\sim 1$. This property
of invariance in the orientational ordering is related to the fixed aspect
ratio of the core, $L/\sigma$, even though the roundness polydispersity is
varied: note that the sine and cosine terms in the excluded area given by Eqn.
(1), are weighted only by $L$ and $\sigma$.

Figure 9: (a) Packing fraction $\eta$ vs. polydispersity $s$ for HRR with core
length and width equal to $L=3$ and $\sigma=0.8$ respectively, and with a mean
disk diameter $l_{0}=1$ (giving mean aspect ratio and roundness
$\kappa=2.\overline{2}$ and $\theta=0.\overline{5}$ respectively. The shadow
and cloud (for both I and N) coexistence curves are correspondingly labeled.
The inset show the fractionation at the I and N shadow coexisting phases
measured through the quantity $\langle l\rangle_{f^{(s)}(l)}/l_{0}$ as a
function of s. (b) Order parameters $Q_{2}$ and $Q_{4}$ as a function of
polydispersity along the cloud and shadow coexisting curves. Inset:
Coexistence scaled densities $\rho^{*}=\rho\langle a\rangle$ as a function of
polydispersity $s$ corresponding to the cloud-shadow equilibrium for both I
and N phases; the labeling is the same as that of panel (a).

Figure 10: Phase diagrams of HRR with $\theta=0.3$ in the plane $\Delta\eta$-$\kappa$, where $\Delta\eta=\eta-\eta^{*}(\kappa)$. $\eta^{*}(\kappa)=a\eta+b$ ($a$, $b$ are constants) is the straight line joining the T-N and I-N tricritical points. The polydispersity coefficients are (a) $s=0$, (b) $s=0.5$, and (c) $s=1$. Different coexistence curves and regions of stability of different phases are labeled. Dotted lines correspond to second-order phase transitions. $\theta$ | 0 | 0.3 | 0.3 | 0.3
---|---|---|---|---
$s$ | 0 | 0 | 0.5 | 1
$\kappa_{t}^{\rm(TN)}$ |  1.940 |  1.676 |  1.669 |  1.644
$\kappa_{ec}$ | 2.210 | 1.861 | 1.857 | 1.844
$\kappa_{t}^{\rm(IN)}$ | 5.440 | 3.915 | 3.932 | 3.976
Table 1: Tabulated values of the mean aspect ratios at the I-N
($\kappa_{t}^{(\rm IN)}$) and T-N ($\kappa_{t}^{(\rm TN)}$) tricritical
points, and also at the critical end-point ($\kappa_{ec}$) for different
values of polydispersity $s$ and fixed value of roundness, $\theta=0.3$. The
corresponding values for HR ($\theta=0$) are included.

Figure 11: Phase diagrams of HRR with $\theta=0.3$ in the plane
$\rho^{*}-\kappa$, where $\rho^{*}=m_{0}\langle a\rangle$ is the scaled number
density. The polydispersity coefficients are (a) $s=0.5$, and (b) $s=1$.
Different coexistence curves and regions of stability of different phases are
labeled. Dotted lines correspond to second-order phase transitions.

Next we describe the changes in the phase behavior of a fluid, whose one-
component counterpart exhibits a first-order I-N phase transition, when
polydispersity is switched on. The particle geometry was chosen such that
length and width are $L=3$ and $\sigma=0.8$, respectively, while the mean disk
diameter is $l_{0}=1$. The mean aspect ratio and roundness result in
$\kappa=2.\overline{2}$ and $\theta=0.\overline{5}$, respectively. Coexistence
packing fractions are plotted as a function of $s$ in Fig. 9(a). In the limit
of zero polydispersity, the shadow and cloud (I or N) phases coincide, as it
should be. Note how the shadow-I and shadow-N curves cross each other at
$s\simeq 0.4$, and also they cross the cloud-N and cloud-I curves,
respectively, at $s\simeq 0.6$. As was pointed out before, this behavior is a
direct consequence of the strong fractionation effect whereby the shadow-I
phase is enriched in particles with larger roundness. It is also clear from
the figure that the coexistence gap is enlarged so that the first-order
transition becomes stronger. In the inset of 9 (a) the scaled mean disk
diameter, averaged with respect to both I and N shadow coexisting phases, is
again plotted as a function of $s$. Clearly, fractionation is strongly
enhanced by polydispersity. The inversion in packing fraction of the
coexisting I and N phases described above is not visible in the number
density, as shown in the inset of panel (b) where the properly scaled
coexistence densities are plotted as a function of polydispersity. Densities
display the usual behavior: the shadow-I (shadow-N) coexistence phase has a
lower (higher) density than the cloud-N (cloud-I) phase. Therefore, the
inversion in packing fraction is not related with a concentration effect but,
as explained above, with the strong fractionation. Finally, panel (b) shows
the coexistence values of the order parameters $Q_{\rm 2n}$ as a function of
polydispersity. They both increase monotonically with $s$. Cloud and shadow
values are very similar, once more a consequence of the invariant
orientational ordering of particles with a fixed core.

All of the above results pertain to the effect of polydispersity on the phase
behavior of HRR when the core dimensions $L$ and $\sigma$ are both fixed. Now
we describe how the whole phase diagram in the plane $\eta-\kappa$ evolves
with polydispersity. Since $\kappa$ is varied, $L$ and/or $\sigma$ will
change. We fixed a moderate mean roundness parameter $\theta=0.3$ and
calculated three different phase diagrams: (i) that corresponding to the one-
component fluid with no polydispersity (already described in Sec. III), (ii)
that with polydispersity coefficient $s=0.5$ and (iii) that with $s=1$.

The results are shown in Fig. 10. Instead of the packing fraction $\eta$, the
difference between the packing fraction and the straight line
$\eta^{*}=a\kappa+b$ connecting the T-N and the I-N tricritical points is used
for better visualizing the cloud and shadow curves, and the coexisting gap. A
first result is that, for high enough polydispersity, the coexistence gap is
enlarged with respect to the one-component case. Also, for $s=0.5$, both
shadow I and N curves are inside the cloud curves, with the packing fraction
of the coexisting I phase (shadow or cloud) being below that of the coexisting
N for any aspect ratio. This is the usual trend. However, for the maximum
polydispersity ($s=1$) the shadow-I and shadow-N curves intersect at
$\kappa\approx 3$, and they also cross their respective cloud curves at
$\kappa\approx 2$ and 2.2, a consequence of the fractionation effect. It is
interesting to note that the T-N coexistence also exhibits strong
fractionation, with the T phase enriched in particles with higher roundness
and lower aspect ratio, while the opposite occurs with the N phase. Note how
the I (shadow or cloud) coexistence curves are always above their N (cloud or
shadow) counterparts in the region of T-N coexistence. When the mean roundness
is moderate, as in the case $\theta=0.3$, the position of the I-N and T-N
tricritical and critical end-points do not change appreciably with
polydispersity, see Table 1. Although the general trend is that polydispersity
enlarges the distance between the tricritical points (the I-N point moving to
the right and the T-N point moving to the right), the difference is visible
only in the second decimal (see table). Plotting phase diagrams in the scaled
density-aspect ratio plane, it is apparent that the shadow and cloud curves
follow the usual trend, with the I or T phases having lower densities compared
with the density of N phase (see Fig. 11). Also, the coexistence gap becomes
considerably larger with polydispersity.

To end this section, Fig. 12 quantifies the fractionation between the
coexisting phases by showing the averaged disk diameter (with respect to the
shadow distribution functions), scaled with $l_{0}$, as a function of $\kappa$
and along the coexistence binodals. As expected, the shadow-I phase has a
large proportion of particles with high roundness as compared to the cloud-N
phase (both cloud phases has this magnitude fixed to 1). The opposite trend is
exhibited by the shadow-N, which is enriched in less rounded particles. This
general trend holds for both polydispersities, $s=0.5$ and $s=1$, except at a
relatively small interval of aspect ratios, $3<\kappa<4$, corresponding to the
case $s=0.5$, see panel (a), for which the I-N phase transition is relatively
weak

Figure 12: Fractionation vs. $\kappa$ measured through the quantity $\langle
l\rangle_{f^{(s)}(l)}/l_{0}$ along the shadow coexistence curves of HRR (see
keybox), with $\theta=0.3$ and polydispersity coefficients (a) $s=0.5$ and (b)
$s=1$.

##  V Crystalline ordering

In the present study we have not taken into account the stability of
nonuniform phases, but they certainly should be present at high enough
densities, at least in the one-component fluid. In Ref. Escobedo  the authors
study the effect of roundness of hard squares on the stability of the T phase
and find that, if the roundness is larger than 0.3 (a value used here to
predict some phase diagrams), the one-component fluid will exhibit a direct
transition from the isotropic to a crystalline phase. We can use a simple
argument to extend this result to HRR and estimate the critical roundness
parameter beyond which a crystal phase is not expected for rounded rectangles
(note that an added polydispersity, not contemplated in this section, would
tend to destabilize the crystal phase regardless of the value of roundness).

To estimate this maximum roundness, $\theta_{\rm max}$, at which the one-
component HRR-fluid with particles of a given aspect ratio, $\kappa>1$,
destabilizes with respect to the appearance of the crystalline ordering, we
use the following procedure. For $\theta_{\rm max}\approx 0.3$, the
destabilization value corresponding to rounded hard squares Escobedo , the
difference between the excluded area (scaled with particle area) of hard
rounded squares in parallel ($\Delta\phi=0$) or perpendicular
($\Delta\phi=\pi/2$) configurations (note that they are identical) and its
maximum value at $\Delta\phi_{\rm max}$ is equal to one half of the scaled
difference but at zero roundness. For rectangles we use the same criterion,
extended to any aspect ratio, and write

$\displaystyle\frac{A_{\rm excl}(\Delta\phi_{\rm max},\kappa,\theta_{\rm
max})-A_{\rm excl}(\pi/2,\kappa,0)}{a(\kappa,\theta_{\rm max})}$
$\displaystyle=\frac{A_{\rm excl}(\Delta\phi_{\rm max},\kappa,0)-A_{\rm
excl}(\pi/2,\kappa,0)}{2a(\kappa,0)}.$ (30)

This equality allows to find the maximum roundness, $\theta_{\rm
max}(\kappa)$, as a function of $\kappa$ for
$1<\kappa\leq\kappa_{c}(\theta_{\rm max})$. We define $\kappa_{c}(\theta_{\rm
max})=\kappa^{*}-(\kappa^{*}-1)\theta_{\rm max}$ as the aspect ratio
corresponding to the intersection between the I-N and T-N spinodals (beyond
which the T phase is no longer stable). Note that, in Eqn. (30) we are using
the excluded area evaluated at the T-like configuration, i.e. for
$\Delta\phi=\pi/2$.

The function $\theta_{\rm max}(\kappa)$ is plotted in Fig. 13. The dotted line
is the straight line $\kappa_{c}=\kappa^{*}-(\kappa^{*}-1)\theta$ for
$0\leq\theta\leq\theta_{\rm max}$. We can see that the maximum roundness,
although slightly higher, has approximately the same value, for any aspect
ratio, as for hard squares. The curves shown in the inset of Fig. 13
correspond to the packing fractions of the I-T bifurcation, evaluated at
$\theta=0$ and $\theta=\theta_{\rm max}$. Also plotted (dotted line) is the
packing fraction at the intersection between I-T and I-N spinodals. The region
between these curves enclose the region of T phase stability as $\theta$
varies from 0 to its maximum value when the crystal phase preempts the T
phase. Of course, this result would be the one resulting from any DFT whose
uniform-density limit gives the SPT. As is well known, this theory
overestimates the packing fraction at which liquid-crystal and nonuniform
phases begin to be stable, especially for small aspect ratios. More
sophisticated theories, with the inclusion of three-body of higher
correlations three-body , are necessary to describe quantitatively the phase
behavior at these aspect ratios. However we are confident that the qualitative
description (except for the precise packing fraction location) of the phase
behavior of HRR described in the present study is the correct one.

Figure 13: The function $\theta_{\rm max}(\kappa)$ (see text) and the straight
line $\kappa_{c}=\kappa^{*}-(\kappa^{*}-1)\theta$ for
$0\leq\theta\leq\theta_{\rm max}$. Both curves approximately enclose the
region where the T phase of HRR is stable. Inset: packing fraction $\eta$ at
the I-T bifurcation, evaluated at $\theta=0$ and $\theta=\theta_{\rm max}$
(continuous curves), and packing fraction at the intersection between the I-T
and I-N spinodals (dotted line), both as a function of aspect ratio $\kappa$.

##  VI Conclusions

In this paper we have defined a new particle model, hard rounded rectangles,
to study the effect of roundness on the stability of the T phase and the
character (first vs. second order) of the phase transitions involved in the
phase behavior of the fluid. The first part of the study is devoted to
characterizing the changes in the phase diagram of the one-component fluid,
while in the second part a continuous polydispersity in the roundness
parameter at fixed core lengths is introduced in order to identify novel
trends in the phase behavior of the polydisperse fluid.

For the one-component fluid we have found that the main effect of roundness is
the destabilization of the T phase: the stability region considerably shrinks
in aspect ratio, and moves to higher densities. Also the I-N and T-N
transitions, which are of first order in the hard-rectangle fluid, transform
into second order or become weaker, the latter scenario occurring if the
aspect ratio lies inside the interval defined by the T-N and I-N tricritical
points. This interval dramatically shrinks with roundness parameter,
disappearing altogether in the limit of discorectangles. In addition, we
believe that the crystal phase should not interfere with the above scenario as
our estimations lead to a wide range of values for the roundness parameter
where the T phase should be stable against the crystal.

When polydispersity is added, and if its magnitude and the mean roundness are
large enough, the I-N transition for certain mean aspect ratios changes from
second to first order, with the presence of fractionation in roundness between
the coexisting phases. We quantified this fractionation by measuring the mean
disk diameter at both coexisting phases, resulting in a N phase enriched in
particles with low roundness. Another important effect when polydispersity is
large enough is a packing-fraction inversion: when the I-N or T-N transitions
are of first order, the I and T phases have higher packing fraction than the N
phase. This phenomenon is related to fractionation: due to the enrichment of
the I or T phases in species with high roundness when they coexist with the N
phase, the packing fraction, being a function of the first and second moments
of the disk-diameter distribution function, will have a larger value.

Finally, as regards how the crystalline phase could modify the results
presented in this study, we expect that the inclusion of a large enough
polydispersity will destabilize the crystal phase with the effect of
increasing the threshold value in mean roundness beyond which the crystal
phase becomes stable.

##  Appendix A Explicit equations for coexistence and bifurcation calculations

The set of coexistence equations (20) can be simplified using the following
properties of the function $A_{\rm spt}(\Delta\phi,l,l^{\prime})$ (see Eqs.
(1) and (3)): (i) Sine an cosine functions in the relative angle $\Delta\phi$
are decoupled from the terms that contain $l$ and $l^{\prime}$. This in turn
implies that the double average of $A_{\rm spt}$ with respect to
$\rho(l,\phi)$ is equivalent to the average of the sine and cosine functions
with respect to $m_{0}(\phi)$ (see Eqn. (7)). Using the Fourier expansion (5),
this in turn gives a result which only depends on the Fourier coefficients
$\\{m_{0}^{(k)}\\}$ of the zeroth moment, while the first moment enters only
through its integrated value $m_{1}$. (ii) Moreover, the first functional
derivative of $\langle\langle A_{\rm
spt}(\Delta\phi,l,l^{\prime})\rangle\rangle_{\rho(l,\phi)}$ with respect to
$\rho(l,\phi)$, which it is needed to calculate the function $c_{1}(l,\phi)$,
will also depend on $\\{m_{0}^{(k)}\\}$ (see Eqn. (18)), which are decoupled
from the polydisperse variable $l$. (iii) Finally, the dependence of
$c_{1}(l,\phi)$ on $l$ is a linear polynomial. This is very important because
taking into account the expression for $f(l)$ from Eq. (21), the integration
over $l$ in the coexistence set of equations (20) can be computed
analytically. Defining the new variables $y^{(\alpha)}_{i}\equiv
m_{i}^{(\alpha)}/(1-\eta_{\alpha})$ with $i=0,1,2$ and $\alpha=c,s$ we obtain,
from (20) and all the above properties, the result:

$\displaystyle
y_{i}^{(s)}=y_{0}^{(c)}e^{-\Delta_{1}(\\{y_{j}\\})}\chi_{i}(\\{y_{j}\\})\frac{\displaystyle\int_{0}^{2\pi}d\phi
e^{-\tilde{c}_{1}^{(s)}\left(\\{y_{0}^{(j)}\\},\phi\right)}}{\displaystyle\int_{0}^{2\pi}d\phi
e^{-\tilde{c}_{1}^{(c)}\left(\\{y_{0}^{(j)}\\},\phi\right)}},$ $\displaystyle
i=0,1,2$ (31) $\displaystyle
y_{0}^{(k,s)}=2y_{0}^{(c)}e^{-\Delta_{1}(\\{y_{j}\\})}\chi_{0}(\\{y_{j}\\})$
$\displaystyle\times\frac{\displaystyle\int_{0}^{2\pi}d\phi\cos(2k\phi)e^{-\tilde{c}_{1}^{(s)}\left(\\{y_{0}^{(j)}\\},\phi\right)}}{\displaystyle\int_{0}^{2\pi}d\phi
e^{-\tilde{c}_{1}^{(c)}\left(\\{y_{0}^{(j)}\\},\phi\right)}},$ (32)
$\displaystyle
y_{0}^{(k,c)}=2y_{0}^{(c)}\frac{\displaystyle\int_{0}^{2\pi}d\phi\cos(2k\phi)e^{-\tilde{c}_{1}^{(c)}\left(\\{y_{0}^{(j)}\\},\phi\right)}}{\displaystyle\int_{0}^{2\pi}d\phi
e^{-\tilde{c}_{1}^{(c)}\left(\\{y_{0}^{(j)}\\},\phi\right)}},$ (33)

while, for the cloud phase, we have
$\displaystyle{y_{i}^{(c)}=\frac{\rho_{0}\langle
l^{i}\rangle_{f(l)}}{(1-\eta_{c})}}$. In the preceding equations we have
defined

$\displaystyle\Delta_{0}(\\{y_{i}\\})=\left(L+\sigma\right)\left(y_{0}^{(s)}-y_{0}^{(c)}\right)+\frac{\pi}{2}\left(y_{1}^{(s)}-y_{1}^{(c)}\right),$
(34)
$\displaystyle\Delta_{1}(\\{y_{i}\\})=\frac{2}{\pi}\left(L+\sigma\right)\Delta_{0}(\\{y_{i}\\}),$
(35) $\displaystyle\chi_{i}(\\{y_{j}\\})=\frac{\langle
l^{i}\rangle_{f(l)}}{\left(1+\Delta_{0}(\\{y_{j}\\})l_{0}s^{2}\right)^{i+s^{-2}}},$
(36)
$\displaystyle\tilde{c}_{1}^{(\alpha)}(\\{y_{0}^{(j)}\\})=-\frac{2}{\pi}\sum_{k\geq
1}g_{k}y_{0}^{(k,\alpha)}\cos(2k\phi),$ (37)

and the packing fractions of the cloud and shadow phases can be calculated
from $\\{y_{i}^{(\alpha)}\\}$ as

$\displaystyle\eta_{\alpha}=1-\frac{1}{1+L\sigma
y_{0}^{(\alpha)}+\left(L+\sigma\right)y_{1}^{(\alpha)}+\pi
y_{2}^{(\alpha)}/4}.$ (38)

Note that, in the definition of $\tilde{c}_{1}^{(\alpha)}(\\{y_{0}^{(j)}\\})$,
a term proportional to the pressure $p^{(\alpha)}$ does not appear (as in Eqn.
(16)) because it should be in the numerator and denominator of Eqns.
(31)–(33). As both pressures should be equal at coexistence they cancel. The
pressure within the new variables is

$\displaystyle\beta
p_{\alpha}=y_{0}^{(\alpha)}+\frac{1}{\pi}\left(\left(L+\sigma\right)y_{0}^{(\alpha)}+\frac{\pi}{2}y_{1}^{(\alpha)}\right)^{2}$
$\displaystyle-\frac{1}{2\pi}\sum_{k\geq
1}g_{k}\left(y_{0}^{(k,\alpha)}\right)^{2}.$ (39)

Looking at Eqns. (31) we can see that the variables $\\{y_{0}^{(s)},\
y_{1}^{(s)},\ y_{2}^{(s)}\\}$ are not independent. Dividing Eqns. for
$y_{i}^{(s)}$ ($i=1,2$) by that for $y_{0}^{(s)}$, we obtain

$\displaystyle\frac{y_{i}^{(s)}}{y_{0}^{(s)}}=\frac{1+(i-1)s^{2}}{s^{2i}\left(T\left(\\{y_{0}^{(\alpha)}\\}\right)+\pi
y_{1}^{(s)}/2\right)^{i}},$ (40) $\displaystyle
T\left(\\{y_{0}^{(\alpha)}\\}\right)=\frac{1}{l_{0}s^{2}}+\left(L+\sigma\right)y_{0}^{(s)}$
$\displaystyle-\left(L+\sigma+\frac{\pi}{2}l_{0}\right)y_{0}^{(c)}.$ (41)

Eqns. (40) can be solved for $y_{i}^{(s)}$ ($i=1,2$) as a function of
$y_{0}^{(s)}$ and $y_{0}^{(c)}$ to find

$\displaystyle
y_{2}^{(s)}=\frac{(1+s^{2})}{y_{0}^{s}}\left(y_{1}^{(s)}\right)^{2},$ (42)
$\displaystyle
y_{1}^{(s)}=\frac{1}{\pi}\left[\sqrt{T\left(\\{y_{0}^{(\alpha)}\\}\right)^{2}+\frac{2\pi
y_{0}^{(c)}}{s^{2}}}-T\left(\\{y_{0}^{(\alpha)}\\}\right)\right]$

Finally the function $\chi_{0}(\\{y_{i}\\})$ can be computed as

$\displaystyle\chi_{0}(\\{y_{i}\\})=\left(\frac{y_{1}^{(s)}}{y_{0}^{(s)}l_{0}}\right)^{s^{-2}},$
(44)

which taking into account (A) is also a function of $\\{y_{0}^{(s)},\
y_{0}^{(c)}\\}$.

Taking into account the preceding discussion, we have a total number of
$2(N_{\rm max}+1)$ independent variables $\\{y_{0}^{(c)},\ y_{0}^{(s)},\
y_{0}^{(k,c)},\ y_{0}^{(k,s)}\\}$, where $N_{\rm max}$ is the total number of
Fourier amplitudes used in the truncated Fourier expansion (5). Thus we need
to solve the single Eqn. (31) for $y_{0}^{(s)}$ ($i=0$) and Eqns. (32) and
(33) for the total number of $2N_{\rm max}$ Fourier amplitudes
$y_{0}^{(k,\alpha)}$ of the cloud ($\alpha=c$) and shadow ($\alpha=s$) phases.
Finally the unknown $y_{0}^{(c)}$ can be computed from the equality of
pressures, $p_{c}=p_{s}$, which guarantees mechanical equilibrium between the
coexisting phases. From this equality and Eqn. (39) we can see that the
variable $y_{0}^{(c)}$ can be written as a function of $y_{0}^{(s)}$ and
$y_{0}^{(k,\alpha)}$ by solving a quadratic equation. Note that the $2(N_{\rm
max}+1)$ variables correspond to the case where both, the cloud and shadow
phases have orientational ordering, for example when the N and T phases
coexist. When one of the coexisting phases, say the cloud phase, is I we need
to solve only $N_{\rm max}+2$ equations because $y_{0}^{(k,c)}=0$ for $1\leq
k\leq N_{\rm max}$. We have solved Eqns. (31)-(33) through a mixed Piccard
iteration method, stopped when a prescribed tolerance criterion is achieved,
$\displaystyle\sum_{\alpha=c,s}\sum_{k=0}^{N_{\rm
max}}\left|y_{0,n+1}^{(k,\alpha)}-y_{0,n}^{(k,\alpha)}\right|\leq 10^{-7}$,
where $n$ label the number of iteration.

When a second order I-(N,T) transition takes place, we can calculate the
corresponding packing fraction as follows: In the close neighbourhood of the
instability of the I phase with respect to the N or T phases we can expand
Eqn. (33) with respect to the small quantity $y_{0}^{(k,c)}$ (with $k=1$ and 2
for N and T symmetries respectively) up to first order. Taking into account
the expression (38) for the packing fraction as a function of $y_{i}^{(c)}$,
and the fact that $y_{i}^{(c)}=y_{0}^{(c)}\langle l^{i}\rangle_{f(l)}$, with
$\langle l\rangle_{f(l)}=l_{0}$ and $\langle
l^{2}\rangle_{f(l)}=l_{0}^{2}(1+s^{2})$, we finally obtain

$\displaystyle\eta_{k}=\left(1+\frac{2g_{k}}{\pi\langle
a\rangle_{f(l)}}\right)^{-1}.\quad k=1,2,$ (45)

where the mean particle area is defined as

$\displaystyle\langle
a\rangle_{f(l)}=L\sigma+(L+\sigma)l_{0}+\frac{\pi}{4}l_{0}^{2}(1+s^{2}).$ (46)

The packing fraction at bifurcation from T to N phase can be obtained by
expanding the exponentials of Eqn. (33) with respect to $y_{0}^{(2n-1,c)}$,
the small odd-Fourier amplitudes, up to first order and evaluating the
resulting integrals at the equilibrium even-Fourier amplitudes
$y_{0}^{(2n,c)}$ (which may be quite large because the T phase can have a high
orientational order). We thus obtain the following set of linear equations,
written in matrix form $B\cdot{\bm{t}}={\bf 0}$ with matrix elements

$\displaystyle b_{kn}=\delta_{kn}-\frac{1}{\pi}g_{2n-1}$
$\displaystyle\times\left(y_{0}^{(2(k+n-1),c)}+y_{0}^{(2|k-n|,c)}(1+\delta_{kn})\right),$
(47)

with $\delta_{kn}$ the Kronecker-delta,
${\bm{t}}=\left(y_{0}^{(1,c)},\dots,y_{0}^{(2m-1,c)}\right)^{T}$ and $m=N_{\rm
max}/2$. This system has a nontrivial solution only if ${\cal
B}(y_{0}^{(c)})=\text{det}\left(B\right)=0$, which allows to find
$y_{0}^{(c)}$ (we should bear in mind the notation
$y_{0}^{(0,c)}=y_{0}^{(c)}$) at bifurcation, and from this, and Eqn. (38), the
value of packing fraction.

The numerical procedure to calculate the two-phase coexistence in the one-
component limit can be obtained from the same Eqs. (31)-(33) by setting $i=0$
and taking into account the limit

$\displaystyle\lim_{s\to
0}\chi_{i}\left(\\{y_{j}\\}\right)=l_{0}^{i}e^{-\Delta_{0}(\\{y_{j}\\})},$
(48)
$\displaystyle\Delta_{0}(\\{y_{j}\\})=\left(L+\sigma+\frac{\pi}{2}l_{0}\right)\left(y_{0}^{(s)}-y_{0}^{(c)}\right),$
(49)

while the packing fraction at the I-(N,T) bifurcation is given by (45) setting
$s=0$.

###### Acknowledgements.

Financial support from Grant No. PGC2018-096606-B-I00 (MCIU/AEI/FEDER,UE) is
acknowledged.

## References

  * (1) V. M. Kaganer, H. Mohwald, P. Dutta, Rev. Mod. Phys. 71, 779 (1999).
  * (2) C. Stefaniu, G. Brezesinski, and H. Mohwald, Adv. Colloid Interface Sci. 208, 197 (2014).
  * (3) E. Barry and Z. Dogic, PNAS 107, 10348 (2010).
  * (4) K. Zhao, R. Bruinsma, and T. G. Mason, PNAS 108, 2684 (2010).
  * (5) K. Zhao, C. Harrison, D. Huse, W. B. Russel, and P. M. Chaikin, Phys. Rev. E 76, 040401 (2007).
  * (6) K. Zhao and T. G. Mason, J. Am. Chem. Soc. 134, 18125 (2012).
  * (7) K. Zhao, R. Bruinsma, and T. G. Mason, Nat. Commun. 3, 801 (2012).
  * (8) Z. Hou, K. Zhao, Y. Zong, and T. G. Mason, Phys. Rev. Mat. 3, 015601 (2019).
  * (9) Z. Hou, Y. Zong, Z. Sun, F. Ye, T. G. Mason, and K. Zhao, Nat. Commun. 11, 2064 (2020).
  * (10) H. Liu, Y. Zong, M. Zu, T. G. Mason, F. Ye, and K. Zhao, J. Colloid Interface Sci. 605, 863 (2022).
  * (11) J. P. R. González and G. Cinacchi, Phys. Rev. E 102, 042903 (2002).
  * (12) J. P. R. González and G. Cinacchi, Phys. Rev. E 104, 054604 (2021).
  * (13) T. Vo and S. C. Glotzer, PNAS 119, e2116414119 (2022).
  * (14) J. A. Anderson, J. A. Antonaglia, J. A. Millan, M. Engel, and S. C. Glotzer, Phys. Rev. X 7, 021001 (2017).
  * (15) H. Schlacken, H.-J. Mogel, and P. Schiller, Mol. Phys. 93, 777 (1998).
  * (16) Y. Martínez-Ratón, E. Velasco, and L. Mederos, J. Chem. Phys. 125, 064903 (2005).
  * (17) Y. Martínez-Ratón, E. Velasco, and L. Mederos, J. Chem. Phys. 125, 014501 (2006).
  * (18) K. W. Wojciechowski and D. Frenkel, Comput. Methods Sci. Technol 10, 235 (2004).
  * (19) A. Donev, J. Burton, F. H. Stillinger, and S. Torquato, Phys. Rev. B 73, 054109 (2006).
  * (20) I. Torres-Díaz, R. S. Hendely, A. Mishra, A. J. Yeh, and M. A. Bevan, Soft Matter 18, 1319 (2022).
  * (21) V. Narayan, N. Menon, and S. Ramaswamy, J. Stat. Mech. P01005 (2006).
  * (22) T. Müller, D. de las Heras, and K. Huang, Phys. Rev. E 91, 062207 (2015)
  * (23) M. González-Pinto, F. Borondo, Y. Martínez-Ratón, and E. Velasco, Soft Matter 13, 2571 (2017).
  * (24) L. Walsh and N. Menon, J. Stat. Mech., 083302 (2016).
  * (25) A. Díaz-De Armas, M. Maza-Cuello, Y. Martínez-Ratón, and E. Velasco, Phys. Rev. Research 2, 033436 (2020).
  * (26) P. A. Monderkamp, R. Wittmann, L. B. G. Cortes, D. G. A. L. Aarts, F. Smallenburg, and H. Löwen, Phys. Rev. Lett. 127, 198001 (2021).
  * (27) R. Wittmann, L. B. G. Cortes, H. Löwen, and D. G. A. L. Aarts, Nat. Commun. 12, 623 (2021).
  * (28) Y. Martínez-Ratón and E. Velasco, Phys. Rev. E 102, 052128 (2020).
  * (29) C. Avendaño and F. A. Escobedo, Soft Matter 8, 4675 (2012).
  * (30) S. Mizani, P. Gurin, R. Aliabadi, H. Salehi, and S. Varga, J. Chem. Phys. 153, 034501 (2020).
  * (31) P. Gurin, S. Varga, and G. Odriozola, Phys. Rev. E 102, 062603 (2020).
  * (32) M. A. Cotter and D. E. Martire, J. Chem. Phys. 52, 1902 (1970).
  * (33) G. Lasher, J. Chem. Phys. 53, 4141 (1970).
  * (34) B. Barboy and W. Gelbart, J. Chem. Phys. 71, 3053 (1979).
  * (35) M. A. Bates and D. Frenkel, J. Chem. Phys. 112, 10034 (2000).
  * (36) M. C. Lagomarsino, M. Dogterom, and M. Dijkstra, J. Chem. Phys. 119, 3535 (2003).
  * (37) A. C. D. van Enter and S. B. Shlosman, Phys. Rev. Lett. 89, 285702 (2002).
  * (38) R. L. C. Vink, Phys. Rev. Lett. 98, 217801 (2007).
  * (39) K. Slyusarenko, D. Constantin, and P. Davidson, J. Chem. Phys. 140, 104904 (2014).
  * (40) M. E. Cates, P. Sollich, and P. B. Warren, Adv. Chem. Phys. 116, 265 (2001).
  * (41) P. Sollich, J. Phys.: Condens. Matters 14, R79 (2002).
  * (42) A. Speranza and P. Sollich, J. Chem. Phys. 117, 5421 (2002).
  * (43) Y. Martínez-Ratón and J. A. Cuesta, J. Chem. Phys. 118, 10164 (2003).