{'fraction_non_alphanumeric': 0.027733527819895208, 'fraction_numerical': 0.008540775003358732, 'mean_word_length': 4.510735060814384, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 3, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Microscopic modeling of contact formation between confined surfaces in solution

We derive a Kinetic Monte Carlo model for studying how contacts form between confined surfaces in an ideal solution. The model incorporates repulsive and attractive surface-surface forces between a periodic (2+1)-dimensional solid-on-solid (SOS) crystal surface and a confining flat surface. The repulsive interaction is derived from the theory of electric double-layers, and the attractive interactions are Van der Waals interactions between particles on the SOS surface and the confining surface. The confinement is induced by a constant external pressure normal to the surfaces which is in mechanical equilibrium with the surface-surface forces. The system is in thermal equilibrium, and particles can deposit to and dissolve from the SOS surface. The size of stable contacts formed between the surfaces in chemical equilibrium show a non-trivial dependency on the external pressure which is phenomenologically similar to the dependency of oscillatory hydration forces on the surface-surface separation. As contacts form we find classical phenomena such as Ostwald ripening, coalescence, and primary and secondary nucleation stages. We find contacts shaped as islands, bands or pits, depending solely on the contact size relative to the system size. We also find the model to behave well out of chemical equilibrium. The model is relevant for understanding processes where the force of crystallization and pressure solution are key mechanisms.


I. INTRODUCTION
Solids brought into contact are ubiqutous and dynamic processes and solid contacts are central to tribology [1,2] and the nature of granular materials in general [3]. Physicists often idealize the contact dynamics and study inert surfaces that deform only mechanically, since stateof-the art surface measurement techniques fail to work in the confined environments where the chemical reactions at the interfaces are important [2], and the vast number of simultaneously occurring chemo-mechanical phenomena that depend on the contact topology and stresses makes modeling difficult [1]. Hence the dynamics of reactive solid contacts, has not received the attention it deserves.
Reactive contact dynamics have important applications in processes such as sintering [4,5], where mineral grains stick together after e.g. compaction without liquification occurring at the grain boundaries, fracture healing/crack sealing [6,7], where material in voids and cracks is rearrange in time such that the aperture decreases without the need of a supersaturated solution, the weathering of rocks and concrete [8][9][10][11][12], which is of fundamental interest to building conservation, in addition to metamorphism, diagenesis and weathering in the Earth's crust [13,14].
Biological applications span from the development of the fracture callus in the reparative stage of bone fracture healing [15] to the initial stages of cell membrane fusion [16].
Stress-induced instabilities in reactive solid-solid boundaries are also responsible for the formation and evolution of stylolites [17,18], and recent experimental [19,20] and theoretical [21][22][23][24] studies of frictional interfaces show that the behavior of the microjunctions (i.e. contact points) between the surfaces is crucial when determining the frictional dynamics.
Popular models of how the grain boundary behaves during pressure solution [25] are growth and dissolution with the presence of a confined thin yet stable liquid film [26], and growth causing stabilizing island-shaped contacts between the surfaces [27]. Recent experiments have shown that growth rims during experiments on the force of crystallization [28], and grains that have undergone pressure solution creep [29], have a structured roughness. This is in disagreement with the liquid film model which predict smooth interfaces [26,28]. We therefore hypothesize that attractive surface-surface interactions are important mechanisms in confined crystallization.
We have previously reported on a model without attractive surface-surface forces, which reproduce known thermodynamics for confined surfaces in solution, as well as the pressure solution and the force of crystallization phenomena [30]. However, this model also predicts a smooth interface. We will therefore in this work address the question whether extending the earlier model by adding a Van der Waals-like interaction between the two surfaces is sufficient to produce a structured roughness.
Questions we want to address are under which conditions we can expect stable contacts to form between the surfaces, how these equilibrate and how they appear once equilibrated, and whether these contacts remain stable in systems where the confining surface is displaced due to the force of crystallization or the crystal surface is dissolved due to pressure solution.
The paper is structured as follows: In Sec. II we introduce all aspects of the model such as the different surfaces and the solution, the allowed transitions and their rates, and the interactions used and how their resulting forces arXiv:2006.02129v1 [cond-mat.soft] 3 Jun 2020 are used to maintain mechanical equilibrium. In Sec. III we present the results of how the system equilibrates, how the equilibrium contacts behave,the contact fluctuations, and finally the out-of-equilibrium properties.The final discussions and conclusions in Sec. IV concludes the paper.

II. MODEL
The model consists of a periodic crystal surface placed in an ideal solution confined vertically by a flat inert surface of the same material at a height h l (t) ∈ R. The latter will from here on be referred to as the confining surface. The crystal surface is modeled using a (2 + 1)dimensional periodic solid-on-solid (SOS) surface.
The SOS condition does not allow for overhangs, hence the surface is described by an array of heights h i ∈ Z (in units of bond lengths l 0 ), where i ∈ [0, L×W ] with L and W being the spatial extents of the system. The top-most particles of the crystal surface, from here on referred to as a surface particles, are the only ones that can take part in transitions. The reactive surface area of the crystal is thus A = L × W . The confining surface is kept inert to reduce the complexity the model, since SOS models are not applicable to systems where two opposing surfaces fluctuate in and out of contact with each other (formed contacts could never break).
The liquid surrounding the crystal has a uniform concentration of solute particles, which limits the model to reaction limited systems. Adding a more realistic description of the liquid is possible, but we will here keep the description as simple as possible. The concentration level may vary in time.
Allowed transitions in the system are dissolution of crystal surface particles into solution and deposition of solute particles to the crystal surface.
Particles interact with other particles through nearest neighbor interactions with bond energy E b . The confining surface is subject to an external force with magnitude F 0 in the direction normal to the confining surface. A repulsive force F λ is generated between the surfaces which increases as the separation decreases. This far the model is identical to the one used in our earlier work on the effect of normal stress on confined crystals [30]. Here we will include additional short-range attractive forces f b (i) between the surface particles and the confining surface. We enforce the mechanical equilibrium of the confining surface, that is, the repulsive force always balances the external-and attractive forces. This involves repositioning the confining surface height h l (t). The acting forces are illustrated together with the surfaces in Fig. 1.
For an initial volume V (0) and concentration c(0), the effective number of solute particles is N s (0) = c(0)V (0). Since the system is periodic particles have no means of escaping or entering the system, which means that the total number of particles is conserved. We can therefore keep track of N s (t) by counting the number of dissolved and deposited particles and add it to the initial value. Hence the concentration at a time t is c(t) = N s (t)/V (t).

A. State Transitions and Rates
A deposition can occur at any site given that the confining surface does not block it. A dissolution can occur at any site given that there is an available neighboring site. These restrictions on the deposition and dissolution reactions ensures that the surfaces do not penetrate into one another.
If the neighboring dissolution site is in the solution, the surface particle dissolves and is removed from the surface. If on the other hand the neighboring site is at the crystal surface, the particle slides one lattice length horizontally. This can be interpreted as an immediate dissolution-deposition chain, and is important to include in order to avoid surface particles in regions where the surfaces are in contact becoming static. Horizontal sliding is the only surface-surface transition we allow, since including transitions up or down kink sites in a SOS model is known to cause an anisotropy between vertical and horizontal diffusion [31,32].
We use that the rate of a particle i dissolving is [30] where ν is a frequency factor, ∆G(i) is the free energy gain by removing particle i from the system, k is the Boltzmann constant and T is the temperature. The deposition rate is proportional to the current concentration c as follows: where we for simplicity have used the same frequency factor such that it can be used to set the time scale for the simulations. The system is in chemical equilibrium when c has an equilibrium value c eq at which no net growth occurs. Using a different frequency ν + would cause the system to equilibrate at a different concentrationc eq = c eq ν/ν + .

B. Free Energies
We model ∆G(i) from Eq. (1) as three terms representing three interactions as follows: where the first term is the particle-particle interaction represented by the breaking of n i nearest neighbor bonds with energy E b , the second term represents the particlesurface attraction, and the last term represents the surface-surface repulsion. The arrows indicate forces acting on the confining surface, which are a constant external force F0, a surface-surface repulsive force F λ , and attractive particle-particle Van der Waals-like forces f b . The confining surface height h l (t) is set such that these forces are in equilibrium, or the surfaces are resting on one another. Between the surfaces there is an uniform ideal solution at a given concentration level.
We set the following criteria for the free energies, which determines the shape and interplay of the different interactions: 1. When in perfect contact (no liquid present), the two surfaces produce an infinite bulk structure consisting only of nearest neighbor bonds.
2. The surfaces are made up of the same material.
3. The attractive particle-particle interactions has a Van der Waals-like decay as 1/d 6 [16], where d ≥ 1 is the center-center separation (in lattice units).
The free energy due to the external force F 0 is F 0 h l , where h l is the height of the confining surface. This is analogous to the gravitational potential. This free energy does not directly depend on the crystal surface structure, hence there is no change in free energy associated with F 0 present in the rate calculations.

Surface-surface repulsion
We use the same expression for the repulsive interaction as in Ref. [30] with the exception that it goes to 0 for no separation: where d i ≥ 1 is the center-center distance between the outer particles of the confining surface and surface particle i in lattice units. The free parameter σ 0 represents the ratio between the strength of the surface-surface repulsion and the particle-particle attraction E b . The decay strength λ D is analogous to the Debye length. We will use λ D = 5 throughout this work and vary the value of σ 0 to study the effect of an increasing/decreasing repulsion.
The reason why we set the interaction to 0 for d i = 1 is that we assume that the repulsion originates from the existence of an electrolyte, and if there is no room for liquid, the only interaction is single bonds with energy E b . Moreover, keeping the repulsion active for d i = 1 would violate the first criteria listed earlier concerning the convergence to a bulk material. The decay to 0 repulsion is probably smooth, and the abrupt cut used here is a simplification. This simplification does not cause instabilities since if d i = 1 at any point, we set the total force to 0 (balanced by a normal force). Hence forces at d i = 1 does not need to be evaluated explicitly.
The total free energy due to this interaction is where ζ ≡ 1 − exp(−1/λ D ).

Particle-particle attraction
The surface particles are confined to a lattice and thus always have a bond length separation and consequently share an energy E b if they are nearest neighbors, and 0 else. The position of the confining surface, however, is not restricted to integer bond lengths, hence intermediate separations d i ∈ (1, 2), that is, separations between the Fig. 2. Illustrations of the three cases in Eq. (6) concerning the free energy between a particle on the crystal surface and a particle in the confining surface separated by a distance di in lattice units. The main idea is that this interaction behaves like the nearest neighbor interaction for integer separations. This is illustrated by di = 2 (left) and di = 1 (right) having no bond and a nearest neighbor bond, respectively. The interaction at intermediate separations (middle), indicated by a stippled bond, is not described by the nearest neighbor model. We model this interaction as a 1/d 6 i Van der Waals-like decay, with added corrections to ensure a continuous transition to 0 for di = 2 for both the function and its derivative. The interaction is given by Eq. (7).
nearest and next nearest neighbor, have some value in G WV (i) ∈ (0, E b ) which we assume decays as 1/d 6 i (Van der Waals attraction).
Summarized we may writẽ This attraction is illustrated for all three cases in Fig. (2). We cannot choose G WV (i) = −E b /d 6 i directly, since this would render the free energy discontinuous at d i = 2. We must therefore correct the potential with a shift (as is often done with the Lennard-Jones interaction [33]). A constant shift, however, would still leave the derivative of the free energy discontinuous at d i = 2, hence we introduce a shift in the derivative as well (a linear term in the original expression). The expression now reads which when inserting d i = 1 and d i = 2 produce the correct values E b and 0 such that Eq. (6) is continuous. In Fig. 3 we have compared this expression to the 1/d 6 i function it was designed to resemble, and it is evident that this is in fact the case. The derivative of Eq. (7) with respect to d i is Fig. 3. Comparison between the attractive interaction used in this workGWV(di), and the Van der Waals-like decaying function −1/d 6 i it is designed to resemble, where di is the surface-surface separation at site i in lattice units. It is evident that there is indeed a strong resemblance. The reason for the difference is that the attractive interaction, as well as its derivative, is required to be 0 at a separation of di = 2 (stippled vertical line). This requirement is a consequence of our assumption that this attractive interaction should act identical to the nearest neighbor interaction used on the crystal surface for integer separations. This is also the reason why for separations di = 1 we have a free energy equal to the bond energy E b .
which equals 0 when d i = 2, which means that the derivative of Eq. (6) is continuous as required.
When we remove a surface particle i with a separation d i ∈ [1, 2), that is, the particle has a non-zero attractive interaction with the confining surface, we are guaranteed that the new surface particle i has d i ∈ [2, 3) and consequently no attraction. The change in free energy used in Eq. (3) is therefore The forces on the confining surface is balanced at all times. This gives the equation where F λ ≥ 0 is the generated repulsive force, f b (j) ≤ 0 is the attractive interaction at surface site j, and F 0 > 0 is the magnitude of the external force. These are illustrated in Fig. 1.
If the two surfaces are resting on one another, i.e., d i = 1 for any i, the net force is always 0 due to a normal force associated with the contact. This is a necessary criterion since we do not model the repulsive forces for d i < 1, Fig. 4. A visualization of the separation (d) dependency of the different free energy sources available in the model which affects the confining surface (top panel) and the resulting total force on the confining surface (bottom panel) for a perfectly flat crystal surface using a repulsion strength σ0 = 1 relative to the bond energy E b , and an external force with magnitude F0/E b A = 0.5. The potential from the constant external force becomes F0d similar to that of a gravitational potential, the attractive interactionGWV is shown in more detail in Fig. 3, and the free energy due to repulsion is given by Eq. (5). We observe multiple equilibria: one where the surfaces are resting on each other (d = 1), one peak in the free energy at a short separation where a strong attraction is in equilibrium with a strong repulsion, and one valley where a weaker (or constant) attraction is in equilibrium with a weaker repulsion. We do not include transitions to the unstable peak equilibrium, since we allow only transitions between stable states. The challenge is thus obtaining the second stable equilibrium here located at h λ ∼ 4.96 indicated by the vertical stippled line. The force peak will appear smoother for rough surfaces.
which could have been included by adding a 1/d 12 i term, i.e. use the Lennard-Jones 12-6 potential [33] instead of a pure Van der Waals attraction.
The height at which the confining surface is resting on the crystal surface, h c , is related to the maximum surface height as h c = max(h) + 1, where h denotes the array of all surface heights h i . is the unit less pressure generated by the external force of magnitude F0, σ0 is the magnitude of the repulsive interaction relative to the binding energy E b (attractive interaction). We observe that as we make the repulsion weaker (σ0 smaller), the equilibria disappear since the attractive interaction overcomes the repulsive interaction before it comes out of range at d = 2. In these cases the only valid equilibrium is when the two surfaces are resting on each other, i.e. the surfaces snap into contact due to an overwhelming attraction. A flat surface under these conditions will always be in contact, however, an arbitrary rough surface can produce many different force profiles under the same conditions. Nevertheless, these examples give an important insight into how we can expect the force profile to look, and what it means for the different equilibria.
The force is calculated as F = −∂G/∂h l . The calculation of the repulsive force is straight forward: where we have set the force in contact equal to zero, which is necessary since the potential is discontinuous due to our abrupt cutting of the interaction potential.
Writing the separation as d i = h l − h i , it is clear that ∂d i /∂h l = 1, such that differentiating with respect to d i or h l yields the same result. Using this, the attractive forces become where the derivative is given in Eq. (8). When d i = 1 for any i, i.e. the surface is resting on one another, we assume that the total force is zero. Hence we do not have to be concerned about instabilities at d i = 1. Figure 4 shows the total free energy together with its individual contributions (top panel) and the resulting total force (bottom panel) for a flat surface. From the bottom panel we see two candidates for out-of-contact equilibria (F tot (d) = 0), however, looking at the top panel we see that the one closest to the crystal surface is unstable. In KMC we allow only transitions to stable states, hence we only consider the far equilibrium point h λ and the contact point at d = 1.
Maintaining mechanical equilibrium thus boils down to calculating h λ and deciding whether to choose it or the contact height h c . We split this into two cases: when h l = h c , i.e. when the surfaces are resting on one another, and when they are not.
When the surfaces are not resting on one another, the confining surface is moved to the closest equilibrium height (h c or h λ ) in the direction of the current total force. Given that the position of the confining surface is h l , and the maximum of the force (the peak in Fig. 4) is located at h m , then this rule translates into the following conditions: Condition 2(a) represents the case where the repulsive interaction is too weak to withstand the applied force. In Fig. 5 the total force for commonly used values in this paper is shown, and we see that condition 2(a) occur for high external forces F 0 and/or a weak repulsive interaction strength σ 0 .
If h λ > h c +2, we are outside the cutoff in the attractive interactions, and h λ has an analytical solution [30]: In practice, we calculate this value, and if it is larger than h c + 1, we know it is a valid solution.
If the surfaces are resting on one another, i.e. h l = h c , the net force is 0 by assumption, and is no longer responsible for the dynamics of the confining surface. If the contacts are unstable, dissolution is the primary mechanism for separating the surfaces, and if the contacts are stable and the system is in chemical equilibrium, the attraction is so strong that condition 2(a) applies almost exclusively. However, if the solution is supersaturated, the repulsive energy in the system will steadily increase as particles deposit, and since the repulsion has a longer range than the attractive interaction, the two equilibra may coexist even when h l = h c . In this case we need a criterion that determines when contact bonds break and the confining surface moves from h c to h λ .
This should be expressed in terms of a rate which depends on the energy barrier between the two states, however, we are unable to do this since we here use rates of the form of Eq. (1) which assumes single particle transitions. In other words, if we want an implicit condition for separating the surfaces, we should use the Eyring rate equation [34] directly.
Since this condition does not impact the equilibrium simulations, and the most important part of the separation mechanism is that condition 2(a) stops applying due to a buildup of repulsive energy, the actual condition is not as important as it might seem. We therefore choose to model it based on simple thermodynamical considerations combined with an attempt frequency.
The probability that the surfaces separate in a given attempt we model as a Boltzmann weight using the change in free energy per area between the contact and the separated states where Z = 1 + W b since the weight associated with staying W stay = 1. When we calculate ∆G tot , we do not include the repulsive free energy at h λ for the points in contact at h c . This represents a sort of retardation time needed to form the electric double layers between the newly formed surface areas.
In this paper we fix the attempt frequency to once every A cycles, which we implemented by scaling P b by 1/A using a frequency of 1. This choice ensures that the probability that the surfaces separate in a given time interval will not have an unphysical dependency on the system size, since the fact that we move only one particle per cycle in KMC makes the time step ∆τ KMC ∝ 1/A [35].

III. RESULTS
We will focus on understanding the equilibrium properties of systems in which stable contacts form between the surfaces. We want to understand how these contact form from an initial separated state, which shapes they possess and why, and how they fluctuate in time.
In order to achieve this we first need a proper definition of a contact. Macroscopically it suffices to define a contact as a point where the surfaces rests on one another. Microscopically, however, we need a less binary definition, since the attractive interactions promoting surfacesurface contacts do not necessarily require the surfaces to rest on one another to do so. We will therefore define a contact as a point which is within the range of the attractive interaction, that is, if h l − h i < 2, then site i is in contact. The contact density is then where A is the area of the (flat) confining surface and κ(i) is 1 if h l − h i < 2 and 0 else. We have ρ WV ∈ [0, 1], where ρ WV = 0 represents completely separated surfaces, and ρ WV = 1 represents the case where all surface sites are in range of the attraction. For any realistic choices of repulsion strength σ 0 , the latter scenario leads to the surfaces joining together perfectly. We initialize the system in an a priori known equilibrium state of the system without an attractive interaction [30], which is obtained by setting an initial concentration ln c(0) = (F 0 /E b A − 3)E b /kT and the confining surface to the height given by Eq. (13). Note that for high pressures this initial state may have initial contacts. We then use the method described in Section II A to keep a constant effective number of particles, such that the closed system will equilibrate automatically. The initial crystal surface is random with an average height equal to 0.
All simulations are done using a 30 ×30 surface lattice. We have done sample simulations using a 50 × 50 surface lattice and observed no noticeable change in the results. The mechanical equilibrium calculations are also very CPU-intensive, meaning that if we want to do a thorough analysis of a vast parameter space, keeping the size as small as possible is very favorable. Periodic boundaries are also very forgiving on the system size, so if we were to open up a boundary, caution should be taken to ensure that the boundary effects are still negligible.

A. Equilibration
In this section we investigate how the contacts in the system form as the system evolves from the initial state to its equilibrium state.
In Fig. 6, we present a high and low pressure simulation of the same system. For the low pressure case we observe three stages. In the primary stage the number of contacts increase very slowly, since most of the contacts which form are unstable due to their small size and dissolve quickly. In the secondary stage growth is rapid since here one or more contacts are above the critical size needed to stay stable. Finally we enter a stage where all the contacts have coalesced into a single contact where the average number of formed contacts equals the average number of broken contacts. This process has a clear similarity to kinetic limited growth theory which passes through the same stages towards the stationary state [36].
For even lower pressures the initial surface-surface separation is too large for any stable contact clusters to appear, and the equilibrium state becomes equal to that with no attractive interactions (total separation). For higher pressures initial contacts are much easier to produce since the surfaces on average are closer to one another, and often start out with parts being in contact already. The primary stage is thus skipped.
We also observe that smaller clusters dissolve in favor of the larger ones, which is analogous to Ostwald ripening [37]. This process occurs since the concentration needed to stabilize a contact cluster decreases as the size of the cluster increases because a larger cluster is harder to dissolve. Hence the largest contact cluster has the fastest net growth, which makes the solution undersaturated for smaller contact clusters. Equally sized clusters will, given enough time, coalesce into a single dominant contact.

B. Equilibrium Contacts
In this section we investigate what determines the equilibrium contact size, shape and stability.
In Fig. 7 we show a selection of equilibrium shapes obtained by simulations at different rescaled temperatures kT /E b , rescaled applied pressures F 0 /E b A using σ 0 = 1. We see three classes of shapes occurring: islands for low contact densities, periodic bands for intermediate contact densities, and pits for large contact densities. The shapes appear rougher for higher temperatures, which is expected due to larger surface fluctuations.
It makes sense that the resulting shape should depend on the contact density, since we can imagine an island transforming into a band when it coalesces with its own periodic image, and the pits forming when the bands grow so wide they merge with their own periodic images.
In practice surface dislocations and/or solute diffusion could set a characteristic length scale promoting the existence of separate islands. These islands could merge by forming bands between one another, which in several directions would lead to pit shapes. Nucleation on system boundaries could also lead to half-islands, bands and pits depending on the size of the crystals growing on the boundaries.
It is clear that the equilibrium contact density is a key factor in deciding the behavior of the equilibrium contact. Since the contact density is not a priori known, we need to investigate how the contact density depends on the system parameters. This relationship is shown in Fig. 8.
For each value of σ 0 we observe a domain where ρ WV = 0, i.e. no stable contacts are formed. This domain represents the cases where the surface fluctuations are not large enough relative to the surface-surface Fig. 6. Equilibration from an initial state at time tν = 0 to a state where the contact density ρWV has stabilized for two different levels of applied pressures F0/E b A using σ0 = 1 and E b /kT = 1. The left column shows results at F0/E b A = 0.5, which is considered a low pressure, and the right column shows results at F0/E b A = 1, which is considered a high pressure. The top row shows how the contact density grows with time. For low pressures, we observe a primary stage where the contacts grow slowly, before reaching a secondary stage where new contacts are formed rapidly. This stage eventually ends when we have reached chemical equilibrium. The stages are separated by vertical dashed lines. The second row shows the accumulated number of formed and broken contacts in the system, which clearly reveal the differences between the three stages just described. The final row shows how the contact cluster centeroids move in time. The crosses are vertical projections separated by a constant number of simulation cycles (the time step is not constant). For F0/E b A = 0.5 we see two prominent clusters coalescing into a single cluster as the secondary stage ends. For F0/E b A = 1, the first stage immediately ends, and smaller clusters dissolve in favor of growing one large dominant cluster. separation for stable contacts to form. Since applying more pressure decreases the separation, and increasing the temperature increases the surface fluctuations, we see this domain curving off with higher pressure F 0 /E b A and higher E b /kT . At higher pressures, contact initiation is not limited by surface fluctuations, since the initial separation is low. In this limit, higher fluctuations simply means that the surfaces can fit less tight. Hence for large pressures we observe that increasing E b /kT results in larger contacts. Fig. 7. Equilibrium shapes with a given contact density ρWV at rescaled inverse temperature E b /kT and rescaled external pressure F0/E b A using σ0 = 1. We observe three classes of shapes: islands, bands and pits. The images have been rendered using Ovito [38]. Fig. 8. Plot of the contact density ρWV vs rescaled inverse temperature E b /kT and rescaled external pressure F0/E b A for various repulsion strengths σ0. The latter is relative to the bond energy E b . Each data point is averaged over ten simulations. A white field means no stable contacts were formed. Higher pressures generally means more contacts since the surfaces are closer, and thus we see more contacts at the top of each plot than at the bottom. The reverse occurs when we increase σ0 (more repulsion), which is why the right figure has a lot less contacts than the left. As we heat the system by lowering E b /kT , we observe fewer contacts, which is due to the surface fluctuations becoming very high, which destabilizes contacts. Hence we see a drop in ρWV as we move left in each plot. The surface height fluctuations are those responsible for creating the initial stable contacts, hence cooling the system by increasing E b /kT makes it very hard to initiate contacts and we see a drop in ρWV after a certain point, unless the pressures is so high that contacts are inevitable. If contacts are inevitable, then lowering the fluctuations means that the surfaces can fit more smoothly on top of one another, and we expect an increase in ρWV. Hence we observe a decay in ρWV to the right in each plot for low pressure and an increase in ρWV to the right in each plot for high pressures. For higher σ0, we see lines at certain pressures that decays slower with increasing E b /kT than their surroundings, hinting to a mechanism that favors contacts to form at certain pressures. Fig. 9. Plot of the most frequently occurring equilibrium contact cluster shape for a given rescaled external pressure F0/E b A, rescaled inverse temperature E b /kT and repulsion strength σ0. The latter is relative to the bond energy E b . These data are obtained by pattern recognition of the same surfaces used to produce Fig. 8. It is clear that there is a strong correlation between the contact density ρWV and the stability of the different shapes. The parameter dependency explained in Fig. 8 thus holds for this figure as well, with the exception that some low contact density cases are missing since they do not correspond to any shape.

Contact favoring pressure levels
At high pressures the surface-surface separation is so low that contact clusters always form. An interesting effect which we will investigate in this section is that at certain lower pressures contact clusters appear much more stable than their immediately lower and higher pressure levels. This can be seen as stripes going into the unstable domain (ρ WV = 0) in Fig. 8. A contact favoring pressure level (CFPL) appears because the corresponding far equilibrium point of the confining surface h λ is very close to an integer value. This means that when the crystal surface fluctuates into contact with the confining surface, the binding energy is close to the maximum value E b at a separation of d i = 1, which makes it harder to remove than if the separation was anything else.
For simplicity lets consider a flat crystal surface at an integer heighth. We know from Eq. (13) that the far equilibrium point is if there are no attractive interactions. The condition that this height is an integer value n then translates into the following equation: The CFPLs then become and the distance between two sequential CFPLs becomes These equations are in agreement with the fact that the CFPLs seems further apart for larger σ 0 . In Fig. 10 we show the predicted CFPLs (top row) together with simulations done at a specific value of E b /kT = 1.2 (bottom row). It is clear that there analytical predictions matches very well with the simulations, however, some levels are not predicted. This is expected since we in the derivation completely ignored the attractive interactions.
The reason why the CFPLs appear flat is that if there is a small gap, the surfaces will snap into contact, resulting in the same scenario as if the confining surface was located at an integer value. This also explains why Fig. 10. Plot of the critera in Eq. (17) (top row) and contact density ρWV (bottom row) vs rescaled external pressure F0/E b A for σ0 = 0.5 (left column) and σ0 = 1.5 (right column). The simulations in the bottom row are performed using E b /kT = 1.2. The top row predicts contact favoring pressure levels (CFPLs) when the function takes integer values. The stippled lines are visual aids to identify these levels. We see that Eq. (17) matches very well where the levels are sparse, but when they become crowded we miss out on some. This is expected since when the contact density increases, the attractive interactions, which we completely ignored when deriving Eq. (17), becomes increasingly important.
the contact density appears to possess discrete levels and transitions from one level to another quite abruptly.
If we instead of solving a mechanical equilibrium problem simply fixed the position of the confining surface, we would expect the contact density to oscillate smoothly as the confining surface position was varied. This behavior is phenomenologically similar to oscillatory hydration forces [16] where the energy cost of removing layers of water confined between surfaces oscillates with a period related to the thickness of the layers.

Equilibrium contact shape vs size
The contact density ρ WV does not contain any information about the shape of the contact. However, by comparing the equilibrium densities in Fig. 8 to the most frequently occurring equilibrium shapes in Fig. 9, we see that the two are strongly correlated.
The equilibrium contact shape minimizes the free energy of the contact, and for high values of E b /kT the free energy is dominated by the binding energies.
The free energy due to binding energies is given by the number of available bonds into solution [30], which means we expect square-like contact shapes for a cubic lattice structure. For lower E b /kT we would expect the shapes to fluctuate around these square-like shapes.
For simplicity we will consider only square surfaces with area L 2 and square contact shapes with area A c = L 2 ρ WV . For a square island shape, the relationship between the contact area and the contact perimeter S I is For the bands we have no perimeter associated with the width of the contact, hence the perimeter is simply The pit shapes are inverted versions of the islands. The area of the pit cavity isĀ = L 2 − A c = L 2 (1 − ρ WV ), such that the perimeter becomes In order to predict which of these shapes are the most stable ones, we need to know which of these three cases has the shortest perimeter for a given value of ρ WV . Comparing S I and S B yields which means that for contact densities lower than 1/4, we expect islands. Comparing S B and S P yields which means that for contact densities higher than 3/4, we expect pits.
Equivalently we may say that we expect bands when ρ WV ∈ [1/4, 3/4]. In Fig. 11 we show the number density of an occurring shape as a function of ρ WV , and we see that the limits derived here agree with the simulations. We also see that islands and pits cannot be stable at the same value of ρ WV without bands also being stable, which means that the bands are a transitional state between islands and pits.

C. Contact fluctuations
In this section we investigate how a stable contact fluctuates in time. In e.g. Fig. 6 we observed that the contact clusters were dynamic, and the contact density ρ WV fluctuated around a mean value in time. We measure these fluctuations as the standard deviation of the time series of ρ WV in equilibrium, that is where X t denotes the average of X in time.
In the top and middle row of Fig. 12 these fluctuations and the corresponding value of ρ WV are presented for . For small contacts we get only islands (I), for medium contacts where ρWV ∈ [1/4, 3/4] we get bands (B) as well, and above this limit we get mostly pits (P). These results are obtained by combining Fig. 8 with Fig. 9. The dashed lines are the theoretical stability limits for the bands from Eqs. (23) and (24). It is clear that they agree well with the simulations. various pressures. We see that the fluctuations increase with increasing temperature while the dependency on the pressure decreases. This happens since for high temperature systems, thermal surface fluctuations are dominant. At high pressures the surfaces are always resting on one another, which means that applying more pressure has no effect. Hence the fluctuations stagnate.
Except for high temperature systems (E b /kT = 0.5), we see that larger contact clusters have lower fluctuations. This is to be expected since the stability of a contact depends on its total number of nearest neighbor bonds, and large compact clusters have a low surface to area ratio. Hence these results suggest that the mobility of a contact cluster increases with increasing temperature and decreasing size.  17). The bottom row shows the average contact clusters height ∆Hc /l0 defined in Eq. (26). The crosses represent simulations where no stable contact clusters were formed. Since the thermal fluctuations of the surface increase with temperature, and the contact is part of the surface, the contact fluctuations also follow this trend. We also see that the contact cluster fluctuations decrease as ρWV increases. The jumps in the fluctuations are clearly caused by ρWV transitioning to a different CFPL. Fig. 13. Fluctuations in the number of broken bonds ∆n− and number of gained bonds ∆n+ sampled every 10 000 time steps for various rescaled inverse temperatures E b /kT and rescaled external pressures F0/E b A using σ0 = 1. The vertical dashed lines represent the contact favoring pressure levels (CFPLs) calculated by Eq. (17). As expected, increasing the temperature increases the fluctuations. The fluctuations stagnate at the point where the surfaces are initiated in contact, such that applying more pressure has no effect. We clearly see a correlation between the jumps in fluctuations and the CFPLs.
Fluctuations in the contact cluster size and the surface roughness are correlated since the more the surface heights are fluctuating, the more often we expect the surface to transition in and out of the contact regime. This is the same mechanism that correlated the pressure (initial height) and the contact size ρ WV in Fig. 8, i.e. smaller height fluctuations are required to form contacts if the surface-surface separation is small. This suggests that the height of the contact cluster should have an impact on its fluctuations. Assuming a stable contact cluster exists, we calculate this height as where the first average is over the domain which is in contact (Ω c ), and the second average is over the domain which is out of contact (Ω c ). From the bottom row of Fig. 12 we observe only small changes in H c as the pressure is increased. The reason for this is that even if the initial surface-surface separation is lower at higher pressures, once a stable contact is formed, we are in practice locked into contact for the rest of the simulation. From here on the system dissolves the part of the surface which is out of contact in order to grow the contact cluster to the optimal equilibrium height H c and size ρ WV . These two quantities are balanced since increasing H c increases the energy cost of growing new pillars to increase ρ WV (more particles are needed).
We also observe jumps in the contact fluctuations in Fig. 12. However, looking at Fig. 13, where we plot the fluctuations in the number of gained and broken bonds separately, it is clear that these jumps are caused by the jumps in ρ WV due to transitions to a different contact favoring pressure level.

D. Out-of-equilibrium systems
Here we will give a qualitative description of the typical behavior of the system out of equilibrium.
When the system has converged to equilibrium at some concentration c eq , a supersaturation Ω ≡ c/c eq − 1 can be introduced by setting the concentration to c Ω (t) = c eq (Ω + 1), (27) and keeping it constant regardless of how many particles are dissolved or deposited. Figure 14 shows the four typical behaviors the system possesses out of equilibrium. For dissolving systems (left column where Ω < 0) we see that if the pressure is low (top row), the equilibrium contacts simply dissolve and we end up with minor fluctuations and no stable contacts between the surfaces. Dissolution at high pressures (bottom row zoom-in) causes the initial large equilibrium contact to dissolve steadily as well, however, due to the high pressure, the surfaces never decouple completely. This causes ρ WV to spike in value whenever the layer closest to the confining surface has been dissolved completely, since this enables the confining surface to snap down to the next layer.
For growth (middle column where Ω > 0), we see that if the pressure is too high (bottom row), the surfaces simply grow into a perfect contact at ρ WV = 1. This occurs because the energy cost of breaking the contacts are very high due to the combined effect of a large number of contacts and a high pressure. This makes it so that being at rest on the crystal surface is the only mechanical equilibrium for the confining surface [condition 2(a) from Sec. II C always holds].
For growth at low pressures (center of the top row), we see that the that the surfaces repeatedly separate,i.e. ,ρ WV drops to 0, after which it grows rapidly as particles are favored to stick in the newly formed cavity between the existing contact cluster and the recently displaced confining surface. This newly formed contact has the same shape as the equilibrium contact. When the healing of the broken contact has stagnated, the crystal surface starts to rise again, and the cycle repeats. This produces the jagged profile for growth at low pressures shown in the top zoom-in of Fig. 14. Each repeating cycle lifts the confining surface one lattice unit.
This repeating separation occurs because at low pressures the initial contact is small, such that the system can build up enough repulsive energy to enable the mechanical equilibrium at h λ where the confining surface is not resting on the crystal surface. When h λ is enabled, whether the surfaces separate is controlled by the bond breaking criteria from Sec. II C.
The period of this cycle should therefore depend on two things: how quickly h λ can be enabled, which is controlled by the supersaturation, and the rate of breaking the bonds, which does not depend on the supersaturation, but on the attempt frequency (which we have set to every A cycles). A larger pressure makes it harder to enable h λ and to break the bonds, thus it should increase the period. If we increase the attempt frequency, the surfaces would separate sooner after h λ is enabled, and if we decease it, the separation would happen later, until a point where the surfaces are able to merge completely before they are able to separate.

IV. DISCUSSIONS AND CONCLUSIONS
The results suggest that the force of crystallization [26], which has been observed and studied in various early experiments [39,40] as well as new [28], can cause a lifting of the confining surface even when there are stable contact clusters between the surfaces. This means that we could have potentially large and irregular surface variations in such a system simply due to attractive interactions promoting contacts between the surfaces. This fact is in agreement with recent experimental observations [28]. For dissolution at high pres- Prior to this point the system has been equilibrated by fixing the effective number of particles. The left column shows results for a dissolving system (Ω < 0), the middle column shown results for a growing system, and the right column shows zoom-ins of the black squares in the respective row. Dissolution at low pressures is simply a decay of the initial contact cluster, and growth at high pressures simply results in the two surfaces merging completely (ρWV = 1). Growth at low pressures (top zoom) periodically causes the confining surface to separate from the crystal surface. For dissolution at high pressures (bottom zoom) we observe an initial stage where the contact cluster dissolves, after which we observe the contact density to spike every now and then.
sures we found the surfaces to possess a continuous state of contact after the dominant equilibrium contact cluster dissolved. Hence the results also suggest that pressure solution might occur with parts of the surfaces being in contact. This is in agreement with experiments observations [29]. This means that we by including attractive surfacesurface interactions made the original model [30] go from producing solely flat interfaces, to producing structural roughnesses in the limits of low pressured growth and high pressured dissolution. This suggests that our hypothesis stating that attractive surface-surface interactions are essential mechanisms in confined crystallization is feasible.
Generally we observe a dominant equilibrium contact cluster forming between the surfaces from an initial state with no such contact. The pressure and temperature dependency of the contact cluster size indicate that it is limited by thermal surface fluctuations (roughness). If the fluctuations are insufficient then no stable contacts are formed. During equilibration we recognize known concepts from kinetic limited growth theory such as primary and secondary nucleation stages, coalescence and Ostwald ripening. These effects are known to appear in island nucleation and step growth [36,37], and their presence suggests that the kinetics are properly treated in our model.
Since the initial surface-surface separation depends smoothly on the external pressure, and the thermal fluctuations depend smoothly on the temperature, the fact that fluctuations limit the contact size tells us that the coexistence of the repulsive and attractive interactions is stable not only for a selective few parameters, but for all. The fact that certain pressure levels promote the formation of contacts due to the initial surface-surface separation being such that formed contacts are especially hard to dissolve, is phenomenologically similar to oscillatory hydration forces [16]. This shows that the model produces effects known to be associated with confined surfaces.
The size of the contact cluster relative to the system size, which we refer to as the contact density ρ WV , was found to be the key parameter deciding whether the equilibrium contact cluster was shaped as an island, a band or a pit, as well as governing the stability of the contact in time (e.g. the fluctuations and mobility). The stability regions of these shapes were in excellent agreement with theoretical predictions. This demonstrates that the model can be used to study details regarding both when and how contacts form between surfaces.
Possible extensions of the model include introducing a Hamaker constant [16] in the attractive interaction to model a different material in the confining surface, and using the Eyring equation [34] for the rates (i.e. calculate energy barriers), such that a parameter for the bond breaking frequency would not be necessary. Moreover, a Lennard-Jones potential [33] could be used between the surfaces instead of the pure Van der Waals term we have used, which would remove the need to distinguish between the surfaces resting on one another and being separated. Using discrete solute particles would enable the study of diffusion limited systems [41]. Elastic interaction could also be added to the surfaces [42][43][44][45].
Details aside, it is fascinating how much interesting physics came out of simply extending the previous model [30] by counting the confining surface as a neighbor. This leads us to believe that there is something simple yet fundamentally correct with our description of how attractive and repulsive forces work together in a confined system.

ACKNOWLEDGMENTS
This study was supported by the Research Council of Norway through the project "Nanoconfined crystal growth and dissolution" (No. 222386). We acknowledge support from the Norwegian High Performance Computing (NOTUR) network through the grant of machine access.