Adsorption superlattice stabilized by elastic interactions in a soft porous crystal

We numerically show that molecules adsorbed in a soft porous crystal form a superlattice (SL) stabilized by elastic interactions. In a mechanically flexible honeycomb lattice model, when the elastic interactions between the next nearest neighboring lattice sites are strong, a long-range ordered 1/3-filling SL state emerges. By calculating the thermodynamic stability, it is found that the SL state is robust against thermal fluctuation. Our results provide a mechanism of elasticity-driven SL formation, which can be utilized for controlling the distribution of adsorbed molecules.

Another consequence of mechanical flexibility is the controllability of the spatial distribution of guest adsorbates.In particular, some SPCs exhibit superlattice (SL) formation upon molecular adsorption.In IRMOF-74-Vhex [19,20] and Co 2 (dobdc) [21,22], it was shown by X-ray diffraction measurements and molecular dynamics (MD) simulations that heterogeneous lattice distortion couples with heterogeneous guest distribution, resulting in the formation of 2 × 2 SL structure.Although guestinduced local framework distortion has been studied, however, the role of long-range guest-guest interaction mediated by the framework's elasticity remains elusive.As shown in other condensed matter exhibiting SL structure, such as nanoparticle quantum dots [23][24][25][26][27], plasmonic SLs [28], and phononic crystals [29], long-range elastic interaction determines the mesoscopic structural formation [25,[30][31][32][33]. In these substances, controlling SL structure is important in enhancing thermoelectric, optoelectric, and phononic properties [23,27,29].Thus, it should be informative also in SPCs, both from condensed matter physics and practical applications, to reveal the role of elasticity-mediated guest-guest interaction leading to SL structure.
To examine the long-range nature of the elastic interaction in SPCs numerically, vast computational costs are required if guest molecules and a deformable host ma-trix are fully incorporated.Thus, a coarse-grained lattice model is needed to elucidate the role of elasticity in SPCs with manageable computational costs.The present authors have constructed a coarse-grained square lattice model, incorporating the adsorption-induced lattice expansion/contraction and hardening/softening [34].Elucidating the connection of the spatial guest distribution with thermodynamics, it has been found that spatial heterogeneity in the stiffness of host frameworks (elastic heterogeneity) leads to the hysteretic adsorption-desorption transition.Extending the coarse-grained model to the honeycomb lattice, where SL structures are observed experimentally [19,21,22], allow us to investigate the mechanism of SL formations in SPCs.
In this letter, we present that an adsorbate SL structure is stabilised by elasticity-mediated guest-guest interaction in SPCs.We consider a coarse-grained honeycomb lattice model which incorporates the adsorptioninduced lattice expansion and hardening.The lattice sites interact with the nearest neighbor (NN) and next nearest neighbor (NNN) sites via simple spring potential.We reveal that elastic heterogeneity leads to the robust hysteresis as well as the square lattice case [34] but, in contrast, a √ 3× √ 3 (1/3-filling) SL state emerges, whose structure is different from one observed experimentally [19,21,22].The SL state is stabilized by the NNN elastic interactions.When the NNN elastic interaction is sufficiently strong, the adsorption fraction exhibits a 1/3 plateau during the adsorption process in a certain parameter region, while there is no intermediate plateau during the desorption process.Correspondingly, the free-energy landscape against the fraction of adsorbed sites takes a local minimum at 1/3, which implies that the SL state is robust against thermal fluctuations.Our findings provide a physical mechanism to realize the SL structure, which can be utilized for controlling the spatial distribution of adsorbed particles.
We construct a flexible two-dimensional honeycomb lattice, as shown in Fig. 1  of the NNN interaction to the NN interaction.Hereafter, we adopt l 0 and k 0 l 2 0 as units of length and energy, respectively.Thus, the potential energy of the unit cell is given by , where {r i∈ } represents the positions of the lattice sites forming hexagon .Each hexagonal unit cell can accommodate only one guest particle.Then the hexagon favors expanding isotropically by the interaction between the guest particle and host matrix, as shown in Fig. 1 (b).This interaction is expressed as additional potential energy where k is the relative energy scale of the guest-host interaction, and α is a swelling parameter, as shown in Fig. 1 (c).Thus, the equilibrium lattice constant and rigidity of a hexagon adsorbing a guest particle are 1 + kα/(1 + k) and 1 + k, respectively.Each hexagon expands/contracts by adsorbing a guest if kα is positive/negative.The sign of k determines whether hexagons harden or soften.In this study, we use a fixed parameter set k = 3 and α = 0.4.
In this study, we adopt an osmotic ensemble [35], whose control parameters are the temperature T , the chemical potential of the guest particle adsorption µ, the number of the host matrix site N host , and the hydrostatic pressure P .The osmotic grand potential is defined as Ω = U − T S + P V − µN ads , where U is the internal energy, S is the entropy, V is the volume, and N ads is the number of adsorption particles.In this study, we fix N host = 2L 2 , where L is the linear system size of the honeycomb lattice.Thus, the lattice site positions r i (i = 1, 2, ..., 2L 2 ) on the honeycomb lattice with periodic boundary condition and guest variables on the hexagons σ ( = 1, 2, ..., L 2 ) taking 1 (presence) or 0 (absence) follow the Hamiltonian given by, Note that the guest variables are on the dual lattice of the honeycomb lattice, i.e., triangular lattice.Hence, the maximum number of adsorbed particles is half of N host .We perform standard Monte-Carlo (MC) simulations for L =12-72 to investigate the hysteretic adsorptiondesorption transition and multicanonical MC simulations using the Wang-Landau (WL) method for L = 12 to study the equilibrium phase transition and free-energy landscapes (see also Supplemental Material).
Figures 2 (a) and (b) show the µ-T phase diagrams at P = 0.5 with the ratio of the NNN interaction to the NN interaction g = 1 and g = 2, respectively.In the case of g = 1, adsorbed and desorbed states are realized in the red and blue regions, respectively, independent of the simulation protocols.Hysteretic behavior is observed in the yellow region.The equilibrium phase boundary is obtained by the WL simulations, crossing the middle of the hysteretic region.Thus, the adsorbed (desorbed) states are not thermodynamically stable below (above) the equilibrium phase boundary in the yellow region.The phase behavior is quite similar to the square lattice case (Ref.[34] Fig. 2).In contrast, in the case of g = 2 (Fig. 2  The SL state stably holds as long as the system does not deviate from the light and dark gray regions by heating or decreasing µ.Thus, hysteresis between desorbed and SL states is observed.Now, we investigate the reason why the 1/3-filling SL state is stabilized.Because the elastic energy depends on configurations of adsorbed sites, effective interactions between guest particles mediated by the elasticity of the host matrix emerge.Fig. 3 (a) displays the enthalpy increase of each local configuration per site ∆h from the ground state at (µ, P ) = (0, 0.5) for g = 1 (The case of g = 2 is presented in Supplemental Material Fig. S1 (ac)).The enthalpy is increased by 2.75 due to an isolated adsorbed site (i).Forming a dimer (ii) or trimer (v, vi) cluster, the enthalpy increases by 0.1-0.2 per site compared to two isolated adsorbed sites.The trimer tends to be isotropic because the adsorbed cluster is harder than surrounding desorbed sites, which is a well-known feature of elastically heterogeneous systems, called Eshelby's argument [33,36,37].On the other hand, when two (iv) or three (vii) adsorbed sites are located on the NNN sites, the enthalpy decreases by 0.05-0.1 per site compared to two isolated adsorbed sites.This implies that the effective interactions induced by the host's elasticity give rise to the formation of the 1/3-filling SL.In more detail, the contribution of the enthalpy change can be divided into the NN elastic interaction, the NNN elastic interaction, and the global pressure term, P ∆V .As shown in Fig. 3 (b), the NN elastic interaction favors forming a connected cluster (ii, v, vi).However, as shown in Fig. 3 (c), the NNN elastic interactions raise the energy when adsorbed sites form a cluster and lower the energy when they are placed at NNN sites.Thus, the NNN elastic interaction is responsible for the formation of the 1/3-filling SL adsorbed state.We note that local configuration (iii), a base of 2×2 (1/4-filling) SL state observed in experimentally [19,21,22,29], is not favored by either NN or NNN elastic interaction.This tendency does not change in the guest-induced lattice contraction case (see Supplemental Material Fig. S1 (d-f)).
The NNN interaction must be sufficiently large for the 1/3-filling SL state to be globally stable.Unlike the above enthalpy comparison under fixed guest particle number conditions, the number of guest particles in our MC simulations is variable by changing chemical potential µ and temperature T .When µ is small, the desorbed state is the most stable because of its little lattice deformation.On the other hand, as µ is increased, a fully adsorbed state, which has three times the number of adsorbed sites of the SL state, is stabilized.For the SL state to be stabilized instead of these two competing states in the intermediate µ region, the enthalpy difference between the cluster of adsorption sites and the SL alignment must be sufficiently large.Because ∆e NNN is enhanced, the SL state emerges for g = 2, but not for g = 1 (see Fig. 2).
To examine the thermodynamic properties and metastability, we next calculate the temperature dependency of the adsorption fraction =1 σ , which is an order parameter of this system.In the following, we show the MC results for g = 2, corresponding to the phase diagram shown in Fig. 2 (b).Fig. 4 (a) shows the temperature dependency of n ads at µ = 5.0 crossing red, yellow, and blue regions in the phase diagram (Fig. 2 (b)).In the standard MC simulations, the transition between the adsorbed and desorbed states exhibits large hysteresis, quite similar to one observed in the square lattice model [34].In sharp contrast, as shown in Fig. 4 (b), n ads exhibits 1/3 plateau during the cooling protocol at µ = 4.9, 4.8 and 4.7 crossing the dark gray region in the phase diagram.This implies that the 1/3-filling SL state is stable in a certain range of temperature and chemical potential.The metastability of the 1/3-filling SL state can be observed in the osmotic grand potential landscape Ω(n ads ) against the adsorption fraction.The case of µ = 4.9 is shown in Fig. 4 (c).The local minimum shifts from n ads = 0 to n ads = 1/3 as the temperature is lowered and converges below T = 0.25, which agrees with the intermediate plateau in Fig. 4 (b).We note that it is difficult to perform the WL simulation below µ = 4.8 due to the huge numerical cost.It is also observed that a hysteretic behavior between two metastable states, 1/3-filling SL and desorbed states, emerges.Fig. 4 (d) shows the temperature dependency of n ads at µ = 4.3 crossing the light gray region in the phase diagram shown in Fig. 2 (b).One can find that the hysteresis appears in a certain temperature range, from T = 0.142 to 0.176.This implies that there is a thermodynamic free-energy barrier between the 1/3-filling SL state, and the SL state is robust against thermal fluctuation below T = 0.176.Thus, it is confirmed that the 1/3-filling SL state exists in the finite parameter region as a metastable state.
In summary, we have shown that the long-range ordered √ 3 × √ 3 (1/3-filling) SL state emerges when NNN elastic interactions are strong.The elastic energy of the SL arrangement is less than one of the fully adsorbed state, which becomes the dominant effect compared to the chemical potential.We have constructed the flexible honeycomb lattice model, incorporating lattice expansion and hardening, and performed standard and multicanonical MC simulations.The phase behavior of this model has been studied.As a result, the long-range ordered 1/3-filling SL state has been obtained by the cooling simulation.Furthermore, by calculating the thermodynamic stability, it has been found that the 1/3-filling SL is robust against thermal fluctuation.
We discuss the reasons for the difference between the SL pattern observed experimentally [19,21,22,29] and the one obtained by our simulations.In an experimental system, IRMOF-74-V-hex, it has been explained that the lattice distortion is caused by capillary condensation of guest particles in the pores.MD simulations [20] have shown that the lattice distortion forms a structure in which regular hexagonal pores with higher adsorbate density are surrounded by anisotropically distorted pores with lower adsorbate density, which is consistent with the 2 × 2 SL pattern.The surrounding pores form a spiral shape around the regular hexagon.This distortion pattern is favored due to the local mechanical properties of the pore, which depend on how the metal ions and organic linkers are connected.Such an anisotropy is not incorporated into our model.The above argument suggests that the pattern of the SL structure can change depending on the local elastic properties, for example, √ 3 × √ 3 SL is formed if the isotropic expansion is favored with strong NNN elastic interactions, and 2 × 2 SL is formed if anisotropically twisted strains are favored.Finally, we discuss the possible applications of adsorbate SLs.Unlike the SLs in other condensed matter [23][24][25][26][27][28], the number of adsorbed molecules is variable; then the adsorbate SLs in SPCs undergo the transitions to homogeneous states depending on temperature and chemical potential.This property is utilized for the switching of the magnetism and conductivity of MOFs as well as optical and phononic properties.In magnetic MOFs, the super-exchange interaction between magnetic moments is modulated by the adsorption of oxygen or nitrogen molecules.Magnetic switching utilizing this modulation has been proposed both experimentally and theoret-ically [38,39].Thus, the formation of intermediate SL states can be utilized for multi-step magnetic switching.Furthermore, some MOFs exhibit semiconductor behavior [40].The periodic potentials induced from the SL may modulate the conductor band structures of semimetal MOFs, which can lead to the transition to metals or insulators.Thus, the adsorbate SL states are expected to be utilized for magnetic/electric sensor devices.
This work was supported by KAKENHI Grant No. JP20H05619 from MEXT, Japan.
Supplemental Material for "Adsorption superlattice stabilized by elastic interactions in a soft porous crystal"

DETAILS OF MONTE CARLO SIMULATIONS Average swelling ratio and volume
Because the average size of the simulation cell varies with the gas adsorption/desorption, the system volume V also varies during the simulations.Therefore, we impose periodic boundary conditions in directions along primitive translation vectors, a 1 = √ 3ℓ 0 (1, 0) and a 2 = √ 3ℓ 0 (1/2, √ 3/2), such that the system becomes periodic under the translation of √ 3ℓ 0 × aL.Here, a = V /V 0 is the average swelling ratio, and V 0 = 3 √ 3L 2 ℓ 2 0 /2 is the reference system volume.

Unit Monte Carlo step
The unit Monte Carlo (MC) step consists of one Metropolis sweep for the adsorption/desorption of guest particles {σ □ }, L iterations of Metropolis sweeps for the lattice sites {r i }, and L iterations of Metropolis updates for the affine displacement of the system to change a = V /V 0 .Thus, the unit MC step consists of N + L × N + L Metropolis updates.The updates of r i and a are restricted to |∆r i | < 0.1 and |∆a| < 0.01, respectively.During the updates of the lattice sites {r i }, displacements causing bond intersects are prohibited to preserve the honeycomb lattice configuration without folding.

Standard Monte Carlo simulation
To capture the hysteretic behavior of the adsorption-desorption transitions, we perform a standard MC simulation, where the temperature T is varied quasistatically.In the standard MC simulations, we iterate 10 4 MC steps for equilibration, and subsequently, 10 4 MC steps for obtaining a thermal average at each T and µ ads .Subsequently, we incrementally change the temperature ∆T = 0.01 above T = 0.3 and ∆T = 0.002 below T = 0.3.In the cooling simulations, we prepare the initial configurations at a sufficiently high temperature T = 1.5.In the heating simulations, we construct two initial configurations: all the plaquettes are occupied by the guest particle, and a 1/3-filling superlattice structure of adsorbed sites is formed.Both elastic energies are minimized.We perform five independent runs for each protocol and evaluate statistical errors; exceptionally, fifteen independent runs are performed at µ = 4.8 and µ = 4.9 for g = 2.

Multicanonical Monte Carlo simulation
To examine the equilibrium phase transitions, the equilibrium probability distribution of the enthalpy H and the adsorption density n = n ads = N ads /L 2 as a function of β = 1/k B T and ν = µ ads /T must be obtained (we fix the pressure P = 0.5 as described in the main text).For this purpose, we adapt multicanonical MC simulations employing the Wang-Landau (WL) method [41][42][43].In the WL method, we calculate the probability distribution function of enthalpy H. Let us consider the probability P (x) of the microscopic state x.If P (x) is uniform, every microscopic state is sampled with equal probability.However, it is computationally inefficient to calculate the probability distribution function of H by uniform P (x).Alternatively, in the WL method, P (x) is proportional to e −g (H) , where g(H) is a weight function.By performing the preliminary run described below, we obtain g(H) ∼ = S(H) + const., where S(H) is the entropy.Hence, the enthalpy histogram A(H) = {x|H(x)=H} P (x) becomes uniform with respect to the enthalpy

ELASTIC ENERGY OF LOCAL GUEST CONFIGURATION
Here, we show in Fig. S5 the enthalpy increase of local configurations from the ground state and the contributions of the nearest and next-nearest elastic interactions at (µ, P ) = (0, 0.5) for parameters that are not shown in the main text.

FIG. 1 .
FIG. 1. Schematic description of the proposed model.(a) Guest particles (dark blue spheres) are adsorbed into the host honeycomb matrix formed by the host particles (orange), inducing isotropic swelling of the host matrix locally.(b) An adsorbed particle strongly interacts with the host particles.(c) Mathematical representation of (a) and (b).Interaction potential V2 arises from the filling of the guest particle in addition to the original host's potential V1.
(b)), a 1/3-filling SL state forming a triangular lattice emerges by crossing the dark gray region in the cooling protocol.The 1/3-filling SL state has a √ 3 × √ 3 unit cell compared to the one of the original triangular lattice, as shown in Fig. 2 (c).The adjacent sites of an adsorbed site are desorbed, and the NNN sites are adsorbed except for defects due to thermal fluctuation.

FIG. 2 .
FIG. 2. Phase diagram of the model at P = 0.5 with k = 3, α = 0.4, and the ratio of NNN interaction to the NN interaction, (a) g = 1 and (b) g = 2.The solid curve represents the equilibrium phase boundary between the adsorbed and desorbed phases, determined from the specific-heat peaks obtained by the WL method.The transition point at T = 0 is analytically determined by comparing the minimum energy of the desorbed and adsorbed states.The boundaries between different colors are determined from the specific-heat peaks obtained by quasi-equilibrium protocols: T1+ and T2+ represent transition points to the desorbed state and √ 3 × √ 3 (1/3-filling) SL state during a heating process, and T1− and T2− represent transition points to the adsorbed state and 1/3-filling SL state during a cooling process.(c) A snapshot of a 1/3-filling SL state at (µ, T, P ) = (4.8,0.2, 0.5) with k = 3, α = 0.4 and g = 2. Light gray hexagons with circles and dark gray hexagons represent adsorbed and desorbed sites, respectively.Short-time averaging over 1000 MCSs with fixed adsorbate distribution is performed to obtain the average lattice distortion.Red lines represent a unit cell of the SL.