Inverse condensation of adsorbed molecules with two conformations

Conventional gas-liquid phase transitions feature a coexistence line that has a monotonic and positive slope in line with our intuition that cooling always leads to condensation. Here we study the inverse phenomenon, condensation of adsorbed organic molecules into dense domains upon heating. Our considerations are motivated by recent experiments [Aeschlimann et al., Angew. Chem. (2021)], which demonstrate the partial dissolution of an ordered molecular monolayer and the mobilization of molecules upon cooling. We introduce a simple lattice model in which each site can have three states corresponding to unoccupied and two discernible molecular conformations. We investigate this model through Monte Carlo simulations, mean-field theory, and exact results based on the analytical solution of the Ising model in two dimensions. Our results should be broadly applicable to molecules with distinct conformations that have sufficiently different entropies or heat capacities.

Conventional gas-liquid phase transitions feature a coexistence line that has a monotonic and positive slope in line with our intuition that cooling always leads to condensation. Here we study the inverse phenomenon, condensation of adsorbed organic molecules into dense domains upon heating. Our considerations are motivated by recent experiments [Aeschlimann et al., Angew. Chem. (2021)], which demonstrate the partial dissolution of an ordered molecular monolayer and the mobilization of molecules upon cooling. We introduce a simple lattice model in which each site can have three states corresponding to unoccupied and two discernible molecular conformations. We investigate this model through Monte Carlo simulations, mean-field theory, and exact results based on the analytical solution of the Ising model in two dimensions. Our results should be broadly applicable to molecules with distinct conformations that have sufficiently different entropies or heat capacities.

I. INTRODUCTION
How a collection of disordered molecular building blocks autonomously arranges into structured and functional materials is of broad interest, from the synthesis of novel materials to our understanding of living matter [1]. Statistical physics provides a powerful framework through which the emergence of structure is connected to the underlying physics of phase transitions. Examples range from the ordering of lipids into membranes and tubules [2] to the drying-induced self-assembly of nanoparticles on a surface [3], the morphologies of which can be captured in a minimal lattice model [4]. Morphologies of thin films of organic molecules deposited on a substrate [5][6][7] play an important role for organic electronics and photovoltaics [8], in particular the control over complex arrested morphologies due to the intricate interplay of thermodynamic and kinetic factors [9,10].
Here we investigate a generic lattice model of adsorption sites, each of which can be occupied by a molecule. Our study is motivated by recent experimental studies of the structure formation of dimolybdenum tetraacetate, Mo 2 (O 2 CCH 3 ) 4 , on a copper Cu(111) substrate in vacuum [11,12]. At room temperature, molecules at submonolayer coverage arrange into aligned chains. Remarkably, upon cooling the substrate the ordered domains partially dissolve and molecules become mobile again, which is reminiscent of inverse melting [13]. The idea of inverse freezing and melting has a long history starting with Tammann [14] in 1903 but has been discussed mostly in connection with rather exotic systems such as helium above 20 bar [15]. Only rather recently it has been realized that adsorbed organic molecules are good candidates to observe inverse transitions more widely [16]. These systems have in common that the liquid-solid phase boundary has an inflection point at high pressure and finite temperature so that further reducing the temperature leads to a reentrance into the liquid phase. Reentrant phase behavior, inter alia, is also found for complex network fluids [17,18], liquid crystals [19], and, generically, liquid mixtures [20,21]. In contrast, here we consider the phase behavior of a single molecular component.
At thermodynamic equilibrium, coexisting phases in an inhomogeneous system have to have the same chemical potential. Changing the intensive thermodynamic variables, along the coexistence line the change of chemical potentials thus has to be equal, dµ g = dµ l , using gas (g) and liquid (l) as a specific example. From the Gibbs-Duhem relation dµ x = −s x dT + v x dp one immediately obtains the Clausius-Clapeyron relation dp for the slope of the coexistence line, where s x is the entropy and v x is the volume per molecule in the corresponding bulk phase, and T and p are conjugate temperature and pressure, respectively. Since we typically expect more available volume in the disordered gas phase together with a larger translational entropy, the slope is positive. A well-known exception is the solid-liquid transition of water exhibiting a negative slope since the volume of ice is larger. Another possibility is that the entropy of the more ordered phase is larger. Clearly, in this case the loss of translational entropy has to be offset by a gain of, e.g., intramolecular entropy. The minimal scenario thus involves two conformations of a complex organic molecule so that in the ordered (solid/liquid) phase the conformation with the larger entropy is found. For dimolybdenum tetraacetate, scanning tunneling microscopy reveals that ordered molecules are standing while freely diffusing molecules lie flat on the substrate ( Fig. 1), thus restricting intramolecular vibrations. Schupper and Shnerb have introduced a spin model that highlights these minimal ingredients necessary for an inverse order-disorder transition [22,23]. It is based on the Blume-Capel model [24,25] with three states for each spin (conformation 1, unoccupied, conformation 2). This model features a tricritical point separating a line of critical points from a first-order transition line. Introducing a degeneracy that endows occupied sites with a higher internal entropy, the line of critical points develops a) Chemical structure of dimolybdenum tetraacetate and the two tentative conformations on top of a copper Cu(111) substrate: (a) Lying molecules binding with ε1. (b) Standing molecules binding with ε2. In addition, intermolecular Mo-O bonding is assumed to contribute a lateral binding energy −ε b , which is absent for lying molecules. Not shown here but important for the model is that both conformations are supposed to have different entropies.
an inflecting point and there is a reentrant behavior back into the disordered paramagnetic phase as temperature is reduced.
Informed by the experiments [12], here we study an even simpler model that qualitatively reproduces the behavior observed in the experiments. As before we consider three states corresponding to unoccupied sites and the two conformations, but now only conformation 2 interacts laterally (Fig. 1). This results in the Ising model for the interacting conformations, which condense into ordered domains below the critical temperature. The non-interacting lying molecules form a gas of mobile molecules outside these domains. Since molecules can freely convert between both conformations, only the total number of molecules is conserved (the so-called semigrand canonical ensemble).

A. Lattice gas model
We consider N molecules adsorbed onto a twodimensional substrate, the atomistic structure of which defines a regular lattice. The substrate is held at constant temperature T . We assume that each molecule can be found in one of two conformations with numbers N 1 and N 2 obeying N 1 + N 2 = N . Specifically for dimolybdenum tetraacetate, conformation 1 is lying on the substrate while conformation 2 is standing upright, cf.
defining the free energies f α (T ). Here, the sum is over all microstates ξ compatible with the conformation α, H(ξ) is the Hamiltonian assigning each microstate an energy, and k B denotes Boltzmann's constant. A coarse-grained configuration of the system is now determined by the positions of the N molecules and their conformations. We assume that molecules bind laterally with binding energy −ε b when both are in conformation 2 and neglect lateral binding otherwise. We assign each lattice site a variablen i = 1 if site i is occupied by a molecule in conformation 2 (standing) andn i = 0 if the lattice site is unoccupied or occupied by a lying molecule (which do not interact laterally). We thus arrive at the free energy functional with h(T ) ≡ f 2 − f 1 the free energy difference of a single molecule in conformation 2 compared to conformation 1, which plays the role of a temperature-dependent external field. We recognize Eq. (3) as the lattice gas (equivalent to the Ising model [26]) with a fluctuating number N 2 of molecules coupled to a finite "reservoir" with N 1 + N 2 = N held fixed.

B. Mean-field theory
The mean-field free energy F derived from Eq. (3) is composed of two terms, the free energy of two ideal gases (corresponding to the two conformations) and an interaction term. We employ as order parameters representing the fraction 0 x 1 of standing molecules and the coverage 0 φ 1, where K is the number of lattice sites (assuming for simplicity that both conformations occupy the same area). The free energy for N non-interacting molecules is given by the entropy to distribute the molecules among K lattice sites with using Stirling's approximation. For the standing molecules there are K sites with coverage φ 2 ≡ N 2 /K = xφ. The N 1 lying molecules can access only the remaining K 1 = K −N 2 lattice sites. Adding both contributions and after some algebra we obtain The interaction term in Eq. (3) is approximated through where we have replaced the inner sum by the average numbern 2 = zN 2 /K of standing neighbors with coordination number z (z = 4 for the square lattice). Adding this term to the ideal gas contribution Eq. (7), the reduced free energy density now reads Here and in what follows we employ dimensionless quantities through rescaling free energies Since molecules can interconvert between both conformations, the first equilibrium condition is which yields the mean-field equation of state Above the critical temperature, T > T c , the curve φ(x; T ) is monotonic. For a fixed φ there is one composition x and the system remains homogeneous. Below the critical temperature, the system becomes inhomogeneous. At coexistence, both the reduced chemical potential and the reduced pressure need to be equal in all coexisting phases. Eliminating φ = φ(x) through Eq. (11), chemical potential and pressure are functions of composition alone with compact expressions Through numerical root finding we determine the compositions x ± obeying µ(x + ) = µ(x − ) and p(x + ) = p(x − ) simultaneously.
To determine the mean-field critical point, we inspect the derivative of the pressure which has zeros for This implies the critical mean-field temperature T MF c = 1 4 independent of h (as expected for the Ising model). Solving x c φ(x c ) = 1 2 with Eq. (11) yields For the critical point to lie within the accessible range, we need x c 1 2 and thus h c = h(T c ) 1 2 .

C. Monte Carlo simulations
We also perform Monte Carlo simulations based on the lattice free energy Eq. (3). Picking a random molecule, we perform one of two trial moves: attempting to move the molecule to one of the z = 4 neighboring lattice sites with probability 1−p conv or switching to the other conformation with p conv . Trial moves are accepted with probability and rejected otherwise, where ∆F = F new − F old is the change of the free energy Eq. (3). Attempted jumps to occupied sites are always rejected. The only parameter of the Monte Carlo algorithm is the probability p conv , which does not influence averages in thermal equilibrium. Most simulation data has been obtained for p conv = 0.1.
To study direct phase coexistence, we employ an elongated simulation box with K = L x × L y lattice sites setting L y = 40 and L x = 3L y . Due to the line tension, the average interface between dense and dilute domains aligns with the shorter box edge so that translational invariance holds along the y-axis but is broken along the x-axis. This facilitates the sampling of x-dependent averages, from which we extract the bulk properties away from the interfaces.

A. Constant free energy difference
We now discuss the phase diagrams following from the analysis of the previous section. To this end, we require an expression for h(T ), and thus for the constrained free energies of the conformations. For our purposes, we decompose where ε α can be interpreted as binding energies with the substrate and the functions s α describe the different entropies of the two conformations (in units of k B ). The mean-field equation of state [Eq. (11)] is plotted in Fig. 2(a) for constant h = 0.35 (i.e., s = 0) and several temperatures T below the critical temperature T c . Following an isotherm through increasing the global coverage φ, the system is a disordered gas of lying molecules with a small fraction x of standing molecules that is slowly increasing. At some φ − we hit the binodal and now dense domains of standing molecules (large φ + and x + ) start to coexist with dilute regions (small φ − and x − ). Further increasing the global coverage φ, these values remain the same but domains occupy a larger fraction κ of the system according to the lever rule Above φ φ + , the system is in the homogeneous dense ("liquid") phase. Increasing the temperature T , the coexistence values approach each other (along the binodal) and the range of global occupations φ − < φ < φ + for which coexistence is possible shrinks and eventually van-ishes at T c . In Fig. 2(b), we plot the coverage φ ± 2 of standing molecules, which is symmetric. We also plot numerical coverages from Monte Carlo simulations, which agree qualitatively but approach each other more slowly at higher temperatures. The critical temperature in the simulations is that of the Ising model, T MC c = [2z ln(1 + √ 2)] −1 0.142. We find that the numerically determined coverages follow Onsager's result for the coexisting densities [27] (22) (obtained in the thermodynamic limit). The standing molecules thus behave as predicted by the standard Ising model independent of h, which controls the total fraction of standing molecules.
We can improve upon our mean-field results through using that the chemical potential of standing molecules at coexistence reads µ * 2 = f 2 − 1/2 in the thermodynamic limit (zero external field in the corresponding spin representation). For the free energy of the lying molecules, we use with fraction ξ ≡ N 1 /K 1 . At equilibrium, the chemical potentials have to be equal. This can be solved for the fraction ξ = [1 + e −(h−1/2)/T ] −1 as a function of temperature T and free energy difference h per molecule. Using the definition we finally obtain a prediction for the total coverage in both the dilute and dense phase. Employing Eq. (22), this prediction is shown in Fig. 2(c) together with the numerical results. While not perfect, the agreement improves dramatically over the simple mean-field ansatz Eq. (8). Finally, in Fig. 2(d) we show the fractions x ± = φ ± /φ ± 2 of standing molecules in both phases.

B. Constant entropy difference
We now lift the degeneracy of the conformational entropies with temperature-dependent h(T ) = ε − T s, where ε ≡ ε 2 − ε 1 is the difference of substrate binding energies and s ≡ s 2 − s 1 the difference of internal entropy of the two molecular conformations. For ε > 0 conformation 2 binds more weakly to the substrate compared with conformation 1. A difference s > 0 implies that the internal entropy of conformation 2 is larger, i.e., more internal modes (vibrations, etc.) are thermally excited. Figure 3(a) shows the mean-field phase diagram for ε = 0.45 and constant s = 0.4, implying the critical value h c = 0.35. While the binodal shifts up, qualitatively we observe a similar behavior to that shown in Fig. 2 for constant h. In Fig. 3(b), we plot the same phase diagram but now for the intensive variables temperature and pressure. The coexistence region becomes a line, crossing which the coverage and fraction jump discontinuously.
In the next step, we make the temperature dependence more pronounced through increasing both ε and s. coverage thus leads to a reentrance into the homogeneous phase. We can easily rearrange Eq. (24) to obtain the total fraction of standing molecules (conformation 2). In Fig. 4(b) we show that this fraction declines and reaches zero at the non-zero temperature T 0 . The system is always homogeneous below T 0 with all molecules being mobile (conformation 1). Specifically, for φ = 1 2 we find T 0 = (ε − 1 2 )/s. For ε < 1 2 we have conventional coexistence without reentrance (T 0 < 0). For s > 0 and 1 2 ε 1 2 + T c s we observe reentrance, whereby at lower temperatures the energetic gain of lying molecules overcomes the entropy surplus of standing molecules. For even larger ε, the system is always homogeneous with a transition from standing to lying molecules at T 0 > T c .

C. Temperature-dependent entropy
Assuming temperature-independent internal entropies s α might not be appropriate for complex molecules. Close to the critical temperature T c , we expand s α (T ) = s α (T c )+(c α /T c )(T −T c ) with heat capacities c α = c α (T c ) per molecule in conformation α (again measured in units of k B ). The free energy difference per molecule now reads with s the entropy difference at T c and c ≡ c 2 −c 1 the difference in molecular heat capacities. Figure 5(a) demonstrates that the resulting phase diagram again exhibits reentrance into the homogeneous phase for sufficiently large c. The total fraction of molecules in conformation 2 drops as before but then rebounds at lower temperatures do to the decrease of the entropy difference [ Fig. 5 Increasing c even further leads to the vanishing of standing molecules and an intermediate homogeneous phase before again ordering at low temperatures. The temperature range T ± 0 of this intermediate phase can be calculated from the vanishing numerator of Eq. (26). Specifically, for s = 0 and φ = 1 2 we find the simple expression with critical value c * 7.05 above which the system goes through the intermediate homogeneous phase.

D. Condensation or freezing?
Before concluding we briefly comment on a subtle point. Both the freezing into a solid and the condensation into a denser liquid phase can be described as orderdisorder transitions. In contrast to condensation, freezing in three dimensions involves the breaking of translational symmetry into a discrete crystal symmetry. Twodimensional systems can still freeze into a solid but without truly long-range positional order, which is destroyed by thermal fluctuations [28]. Strictly speaking, therefore, no crystal exists in two dimensions. While a molecular monolayer comes very close to two dimensions, the underlying substrate induces an external potential that, even if weak, breaks translational symmetry and thus can reintroduce long-range order [29]. Note that a lattice model is insufficient to capture the spontaneous breaking of translational symmetry since the discrete lattice cannot accommodate a continuous symmetry. Rather than freezing into a periodic crystal, our model thus describes condensation of molecules into dense domains.

IV. CONCLUSIONS
The Ising model provides a rigorous statistical mechanics underpinning for a wide class of disorder-order transitions and serves as a paradigm for reversible structure formation. It captures the competition between shortrange isotropic attraction and translational entropy, describing the condensation of a dense liquid coexisting with a dilute gas. Crossing a first-order transition gives rise to phase ordering kinetics that can be exploited to control and create non-equilibrium structures and morphologies [30].
In contrast to atoms and colloidal particles, molecules have internal degrees of freedom. For molecules with approximately isotropic intermolecular interactions, the Ising model is still an appropriate large-scale description but the "field" may acquire a dependence on temperature due to the intramolecular entropy. An additional competition between intramolecular energy and entropy opens a route to reshape the phase diagram. Here we have studied the case of two molecular conformations, whereby intermolecular interactions are only assumed for one confirmation (here conformation 2). This restriction has allowed us to exploit the analytical solution of the Ising model in two dimensions to accurately construct the phase behavior, which is determined by two temperatures: the critical temperature T c of the Ising model below which domains form (driven by the interactions of conformation 2) and the temperature T 0 below which conformation 2 vanishes. While T c is fixed, tuning T 0 allows to shape the coexistence region and to find parameters for which reentrance into a homogeneous disordered phase is possible.
Such reentrance upon cooling has been observed recently for dimolybdenum tetraacetate and has been attributed to two different molecular conformations [12]: Lying molecules (conformation 1) bind more strongly due to the metal-metal interaction (ε 1 < ε 2 ), which at the same time suppresses internal degrees of freedom (s 1 < s 2 ) due to steric hindrance with the substrate (cf. Fig. 1). To relate the scenario of inverse condensation to the experiments of Ref. 12, we first note that (for the ordered phase to be present) room temperature needs to be below the the critical temperature. This condition implies the lower bound ε b > 25.7 meV/(zT c ) ≈ 45 meV. As reference, the effective binding energy of C 60 molecules has been estimated to be ε b ≈ 235 meV [31] and thus is five times larger. The difference in substrate binding energies between the two conformations of dimolybdenum tetraacetate has been determined toε ≈ 410 meV [11], hence dimensionless ε =ε/(zε b ) of order unity seem reasonable. As an alternative scenario, we have considered a difference of heat capacities between both conformations, which also leads to reentrant behavior for differences exceeding a few k B 's. Our results thus indicate that inverse condensation should indeed be observable for a wide range of adsorbed organic molecules.
Beyond the phase behavior of adsorbed molecules, our results are potentially interesting for condensates of proteins forming membraneless organelles inside the cell [32,33], which are studied intensively at the moment. Reentrant behavior has already been reported [34,35] and could be studied theoretically in on-lattice "stickersand-spacers" models for multivalent proteins [36].