# The role of the diffusion in the predictions of the classical nucleation theory for quasi-real systems differ in dipole moment value

Kajetan Koperwas, Filip Kaśkosz, Frédéric Affouard, Andrzej Grzybowski,
Marian Paluch

<figure><img src="image_1.png" /></figure>

## To cite this version:

Kajetan Koperwas, Filip Kaśkosz, Frédéric Affouard, Andrzej Grzybowski, Marian Paluch. The role of the diffusion in the predictions of the classical nucleation theory for quasi-real systems differ in dipole moment value. Scientific Reports, 2022, 12 (1), pp.9552. 10.1038/s41598-022-13715-y. hal-03715087v3

HAL Id: hal-03715087

https://hal.univ-lille.fr/hal-03715087v3

Submitted on 21 Jul 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

Check for updates

OPEN

# The role of the diffusion in the predictions of the classical nucleation theory for quasi-real systems differ in dipole moment value

Kajetan Koperwas<fnref n="1" />✉, Filip Kaśkosz<fnref n="1" />✉, Frederic Affouard², Andrzej Grzybowski<fnref n="1" /> & Marian Paluch<fnref n="1" />

In this paper, we examine the crystallization tendency for two quasi-real systems, which differ exclusively in the dipole moment's value. The main advantage of the studied system is the fact that despite that their structures are entirely identical, they exhibit different physical properties. Hence, the results obtained for one of the proposed model systems cannot be scaled to reproduce the results for another corresponding system, as it can be done for simple model systems, where structural differences are modeled by the different parameters of the intermolecular interactions. Our results show that both examined systems exhibit similar stability behavior below the melting temperature. This finding is contrary to the predictions of the classical nucleation theory, which suggests a significantly higher crystallization tendency for a more polar system. Our studies indicate that the noted discrepancies are caused by the kinetic aspect of the classical nucleation theory, which overestimates the role of diffusion in the nucleation process.

Although the crystallization process is a commonly known phenomenon, the complete understanding of its nature is still far from being achieved. It is mainly due to its complexity, which finally makes that some systems easily crystallize, whereas others do not exhibit any symptoms of crystallization, even at deep undercooling, and finally form the glass. Thus, the complete understanding of this process, including the determination of the physical factors, which govern its occurrence, seems to be a crucial task for contemporary condensed matter physics<fnref n="1" />⁻⁵.

Consequently, through the last decades, various theoretical and computational approaches to study the crystallization phenomenon have been proposed. The computational experiments mainly focus on the possibility of the precise calculation of the order parameter, which enables, e.g., the estimation of the time scale, at which the ordered phase within the liquid system appears for the first time (the mean first passage time method)⁶⁻⁸. Then, the structure of the formed crystal⁹ and the direct evolution of its size, can be immediately monitored¹⁰. However, from the experimental point of view, the theoretical methods employing macroscopic features of the system are of more practical importance. Therefore, a variety of theoretical descriptions for the crystallization process have been proposed¹¹. Most of them are grounded on the same concept, i.e., the crystal phase starts to spread only if the nuclei of a given (critical) size (and/or shape) are formed¹¹. Among those models, the most widely used (probably due to its remarkable simplicity) is the classical nucleation theory (CNT)⁵,¹²,¹³. According to the CNT, the critical nuclei's stability is determined by the surface and bulk contributions to the free energy. Consequently, discussed concept predicts that the crystallization process consists of two stages—the formation of the nuclei of critical size (nucleation) and their growth (crystal grow). The first step can be estimated by the number of the nuclei created within the given volume during the fixed time, i.e., the nucleation rate *N*, whereas the second one, i.e., crystal growth rate *U*, describes the velocity of the growth of the crystal structure within the liquid. As a consequence, the overall crystallization can proceed only when *N* and *U* are coupled. This simple idea enables explanation of essential experimental observations, i.e., it justifies why some systems can be supercooled up to the glass transitions, whereas others crystallize during cooling, and why some supercooled

<footnote n="1">¹Institute of Physics, University of Silesia in Katowice, 75 Pułku Piechoty 1, 41-500 Chorzów, Poland. ²Université de Lille, CNRS, INRAE, Centrale Lille, UMR 8207-UMET-Unité Matériaux et Transformations, 59000 Lille, France. ✉email: kajetan.koperwas@us.edu.pl; filip.kaskosz@smcebi.edu.pl</footnote>
liquids crystallize during the heating from the glass (it is so called cold crystallization)¹⁴. In the first case, the separation of *N* and *U* plays a key role. The *U* curve is located closer to the melting temperature than *N* curve. Hence, at small supercoolings when both components of the crystallization process are substantially separated, *N* is insufficient to create the stable critical nuclei, which would subsequently growth. On the other hand, at a deeper supercooling, the critical nuclei can be created, but then their growth is suspended by the scarce value of *U*. Finally, the substance does not crystallize. At this point, it must also be noted that a slower cooling rate implies that the system persists at given thermodynamic conditions for a longer time. Therefore, the chance for the creation of (at least one) critical nuclei is higher. Hence, CNT considers the effect of the cooling rate as well. In the second case, when a deeply supercooled liquid is heated, the critical nuclei created at a deep supercooling begin to achieve the temperatures at which *U* exhibits high values. Thus, we can observe the crystallization process, which previously, during cooling, was unable to take place. Contrary, if *N* and *U* curves are close to each other, the optimal temperature range for the crystallization process appears. Then the critical nuclei are formed, and subsequently, they freely grow.

In this paper, on the basis of the two highly similar systems, we challenge the prediction of the CNT. Interestingly, despite the fact that at given supercooling one of studied systems exhibits significantly higher values of the *N* and *U*, the crystallization event for this system is not observed. Our examinations suggest that the observed inconsistency between CNT predictions and computational experiment results is caused by the differences in the molecular mobility between studied systems. Consequently, we show that the crystallization process's kinetic aspect should not be straightforwardly linked with diffusion, as CNT assumes.

The CNT has frequently been used for the theoretical description of experimental and computational experiments through the last decades. However, the computational experiments deserve particular attention because it enables to examine the crystallization tendencies on the most fundamental level of intermolecular interactions. For this purpose, the simple model systems characterized by the well-defined intermolecular potentials can be used. The most frequently studied systems are those in which pairwise intermolecular interactions are described by the Lennard-Jones potential or its approximation valid at short distances, i.e., the soft-sphere potential¹⁵⁻¹⁸. The mentioned choice is justified by the fact that the Lennard-Jones potential can be theoretically derived on the basis of the interactions between permanent and induced dipole moments. Consequently, those simple model systems were used to verification of the CNT¹⁹⁻²¹ as well as also to study the influence of the attractive and repulsive intermolecular interactions on the crystallization tendency²²⁻²⁸. Reported studies suggest that the increase in the repulsion results in the decrease in the nucleation barrier and interfacial free energy²⁹. The other examinations focused on the role of the attraction in the crystallization process, deliver the conclusion on the positive impact of the intermolecular attraction on the reduction of the time needed for the crystallization at given temperature³⁰. Simultaneously the different approach, i.e., the computational studies performed on the hard molecules, revealed that the molecular anisotropy ignored by simple model systems is crucial in determining the phase diagram of the system³¹⁻³⁶. However, it must be mentioned that for hard molecules, the temperature enters the thermodynamics only in a trivial way³⁷. Consequently, the alternative models, which consider the interactions between non-spherical molecules have been developed, e.g., Kihara³⁸ potential, the Gaussian overlap model³⁹, and the Gay-Berne potential⁴⁰, and prove the important role of the structural anisotropy in the thermodynamics and dynamics of the studied systems. However, from the experimental point of view, the most natural is the all atom-atom (or site-site) interactions approach, which unfortunately requires much more computational effort⁴¹. Nevertheless, the all atom-atom approach makes that the structure of complex molecules can be reflected, and therefore, the closing agreement with the experiments may be expected. As a consequence, the structural⁴¹⁻⁴⁴ and dynamical⁴⁵⁻⁴⁹ properties of many model system have been deeply examined concluding that this approach can be successfully applied for slightly non-spherical molecules⁵⁰. Following this result, the very recent study reports that the permanent dipole moment orientation within the anisotropic molecules is of crucial importance for the crystallization process. The two analogical systems, which vary exclusively in the orientations of the dipole moment, exhibit entirely different stability behavior despite that the same isobaric conditions and identical cooling rates are applied⁵¹. Briefly speaking, the perpendicular to the longest molecular axis orientation of the dipole moment favors crystallization. In contrast, the deep supercooling of the system with parallel to the longest molecular axis orientation of the dipole moment is easily achieved. This outcome is not only relevant from the experimental point of view, but it is also important for further computational studies because it emphasizes the practical utility of the model systems tested therein. Only slight modifications of the molecular architecture result in drastically different crystallization tendency. Hence, model systems from Ref.⁵¹, which comprise the so-called quasi-real molecules, seem to be promising candidates to examine the crystallization process and then the predictions of the CNT. At this point, it is also worth justifying that the use of the quasi-real molecules, i.e., the molecules which mimic the real ones but cannot exist in reality, helps to eliminate the uncontrolled effects of various intramolecular factors on the considered process or physical quantity. At this point it is also worth mentioning that in contrast to simple models, the results obtained for one of the quasi-real systems cannot be appropriately scaled to reproduce results registered for another system. Thus study on quasi-real systems provides promising alternative to typical computational experiments.

## Results

Similarly, to our previous examinations of the system *I*, we began the studies of system *II* from the constructing the perfect FCC crystal structure constructed from 2048 molecules and heating it from 10 K up to the temperature which is about 50 K higher than the temperature at which we observe a rapid increase in the volume, see Fig. 1.

On the basis of our recent results for the system *I*, we can state that the rapid increase in the volume from about 0.086 to 0.094nm³, which is observed around *T* = 220K, indicates on the melting of the crystal structure. Therefore, the thermodynamic conditions at which the system *II* is in the liquid phase can be recognized in the
<figure><img src="image_4.png" /><figcaption>Figure 1. The temperature dependences of the volume for two studied model systems registered during heating and cooling are shown. The black lines represent the heating of the determined crystal structure. In the insets, the schemes of the structures of RM are presented.</figcaption></figure>

similar way, see increase in volume from about 0.088 to 0.101nm³ at T = 420K. Next, we cooled the systems up to the starting temperatures. During this process the drop in the volume can be observed. The latter indicates on occurrence of the crystallization process, which takes place at Tcr = 130 and 280K, respectively for the systems I and II.

According to CNT the nucleation rate is expressed as follows¹³

$$N = \rho_{liq}^{4/3} \sqrt{\frac{\gamma}{k_B T}} D \exp\left(-\frac{\Delta W}{k_B T}\right), \quad (1)$$

where ρliq is the number density of the liquid, D is a diffusion constant, kB is the Boltzmann constant, and ΔW = (16/3) γ³ / (ΔGv)² is the nucleation barrier, in which ΔGv denotes the driving force per volume unit (i.e., the difference between Gibbs free energy for liquid and bulk phases) and γ is the interfacial free energy (IFE). The next physical quantity determining the occurrence of the overall crystallization process is the crystal growth rate, role of which can be computed by the following expression

$$U = A_U(T) \cdot f(T) \cdot \left[ 1 - \exp\left( -\frac{\Delta G}{k_B T} \right) \right], \quad (2)$$

where A_U(T) describes the molecular mobility and can be approximated by D · a/λ², in which a is the average width of the crystal lattice spacing (a ≈ ρcr, ρcr denotes the number density of the crystal), λ ≈ ρliq⁻¹/³ is the atomic jump distance, whilst f(T) describes the grow mechanism, which for the normal growth ≈ 1. For all thermodynamic conditions at which during the cooling systems remain in the liquid phase the diffusion constant is determined using the GROMACS software on the basis of the mean square displacement calculated for atoms for long times. In this way both translational and rotational contributions to the molecular motion are considered. Subsequently, the estimated dependences D(T) are approximated by the Vogel-Fulcher-Tammann equation for the needs of further analysis (see Fig. 5a, which we discussed later). To estimate the value of ΔG we use the method proposed by Gutzow⁵² with the redefined integration pathways⁵³, according to which, at isobaric conditions, the driving force for the crystallization takes the following form ΔG(T) = −∫TmΔS(T)dT, where ΔS is the difference in the entropy between the liquid and the solid phases. Taking into account the melting boundary conditions ΔSₘ = ΔHₘ/Tₘ, ΔS can be calculated using obtained directly from simulation-runs values of the enthalpy H(T) and the classical relation between enthalpy and the entropy, T = (∂H/∂S)ₚ. The estimation of the melting temperature has been done using the liquid-solid coexistence method. In this order, we visualized the structure up to which each system crystallizes. Then, we determined the fragments characterized by a high degree of order, which for both systems are characterized by the triclinic shape and consist of molecules placed in corners. Based on the latter, we constructed another crystal structure and equilibrate it at the temperature close to Tcr for both systems. Since we observed that the small defects occur again, we selected the set of 5 × 5 × 5 molecules within which the created crystal structures were highly ordered and those crystal fragments are used for further examination. On their basis, we construct the crystal structures consisted of 2250 molecules, and equilibrate it at the temperatures significantly smaller than Tcr. Subsequently, we heat the systems to confirm that the crystals are stable at higher temperatures. The results are presented in Fig. 1 of the main text. It can be seen that created crystal structures do not tend to melt, although, during the heating process, the tiny and step increase in volume is detected for both systems. Probably, these changes of volume are results of the transformation to the different polymorphic form, which is more stable at higher temperatures. This scenario seems to be supported by the evident visible change in the temperature dependence of volume, which is observed around
100K during the cooling of the system II. However, to confirm this suspicion further researches are required. Nevertheless, it is worth mentioning that created structures are more stable at higher temperatures than the ones resulted from starting FCC configuration. This fact encourages that the crystal structures established by us can be used to determine T<sub>m</sub>. Hence, we constructed the special biphasic simulation box, within which 3456 molecules had been equally divided between the crystal and liquid phases separated by a small gap. Due to the fact that at melting conditions, the solid and liquid phases remain in the thermodynamic equilibrium, T<sub>m</sub> can be determined by the examination of the behavior of the biphasic system. The performed simulations of biphasic box last for 5ns and have been done at range of temperatures differ by 1K. Following the visual examination of the obtained configurations we determine that melting temperatures, which equal T<sub>m</sub> = 150K and T<sub>m</sub> = 326K for the systems I and II, respectively. At this point, we have to comment that in our previous studies<sup>51</sup> we determined the T<sub>m</sub> for the system I, using liquid-solid coexistence method as well, and we obtained that T<sub>m</sub> = 194K. However, in that experiment we did not determine the crystal structure. Instead of that, we employed the structure up to which the liquid had crystalized. Consequently, we probably employed the structure stable at higher temperatures instead of the one, which is the most energetically optimal. However, in this work, we intend to focus on the most fundamental case, i.e., we estimate the crystallization tendency against the desired (and the most energetically optimal) structure. In this context, it is worth noting that the values of T<sub>m</sub> estimated herein are in similar relations to those of T<sub>cr</sub> and also to the temperature at which initial crystals melt, which suggests that the determined structures are mutually appropriate.

However, the most challenging is the determination of the γ. Fortunately, the special computational method for calculation of γ have been proposed. The two main approaches are the cleaving potential method<sup>54–56</sup> and the capillary fluctuation method<sup>57–61</sup>. Despite that both methods are applicable only at the melting conditions, they strongly differ in the way of work. In the cleaving potential method, the biphasic solid–liquid system is transformed into two separate systems (liquid and solid) by means of external potentials. Then, γ is estimated on the basis of the work which is performed by those potentials during the transformation process. However, the precise application of this method is associated with some technical difficulties. It is because the reversibility of the transformation process must be ensured, and therefore, the accurate control on the transformation process is needed<sup>62</sup>. Alternatively, γ can be calculated in the more direct way using the capillary fluctuation method (CFM), which requires only one simulation run, during which any knowledge of the complex process of the interface creation from separated bulk systems is not needed. Instead of that, through the simulation run, the fluctuations of the interface are measured. It enables the estimation of the interface stiffness, which is related to γ. The remarkable advantage of CFM is the fact that it considers the anisotropy of the interface, whereas the cleaving method is recognized as more accurate. Till now, both methods have been applied to calculate γ values for model systems such as hard-spheres<sup>55,63</sup> and Lennard–Jones<sup>25,56,60</sup>. It must be however noted that for the real materials the CFM is more often employed, which is mainly stimulated by the ease of its application. Consequently, using the CFM the γ values have been calculated for metallic compounds<sup>57–59,64</sup>, alloys<sup>65,66</sup>, and a few molecular systems<sup>67</sup> including pharmaceuticals<sup>68–70</sup>. Hence we decided to employ CFM to determine γ for studied herein systems. The use of CFM requires creation of the biphasic box. However, the considered solid–liquid interface must be the quasi-one-dimensional, and therefore the special geometrical conditions of the simulation box have to be ensured, i.e., when interface is perpendicular to the length of the system, L<sub>x</sub>, its thickness must be much smaller than its width, L<sub>z</sub> ≪ L<sub>y</sub>. Then the interface fluctuates only in the one dimension (x). Consequently, we construct the box containing 5000 RLM divided equally to the crystal and liquid phases. It is also worth mentioning that due to boundary conditions the simulation of the biphasic system implies the existence of two interfaces, of which the fluctuations magnitudes are studied. The convenient way used to determine the temporary position of the interface is the calculation of the rotational-invariant order parameter<sup>64,71–75</sup> (RIOP) for geometrical center of the molecules. The RIOP enables the distinction between solid-like and liquid-like molecules, because the solid-like molecules are characterized by the significantly higher values of the order parameter. The example of obtained results is presented in Fig. 2a, where the calculated RIOP for each molecule of system I is plotted as a function of the position of the molecule in the dimension perpendicular to the interface plane (x direction).

As we already mentioned the liquid-like and the solid like molecules can be clearly distinct. The evolution of the RIOP can be described by the following function RIOP(y) = (o<sub>s</sub>+o<sub>l</sub>) / 2 + (o<sub>s</sub>-o<sub>l</sub>) / 2 [tanh(y-h<sub>1</sub>(x)) / (δ<sub>1</sub>) + tanh(y-h<sub>2</sub>(x)) / (δ<sub>2</sub>)], where o<sub>s,l</sub> are the average values of the RIOP in the solid and liquid, δ<sub>1,2</sub> are effective widths of the interfaces, and h<sub>1,2</sub>(x) are functions describing the positions of the interfaces in capillaries, i.e., sections from x to x + Δx, which are orthogonal to the interface. During the simulation run the h(x) describes the interface fluctuations. The latter can be Fourier-transformed leading to the following expression for power spectrum of the quasi-one dimensional interface ⟨|h(q)|<sup>2</sup>⟩ = k<sub>B</sub>T<sub>m</sub> / (L<sub>x</sub>L<sub>z</sub>γ<sub>m</sub>q<sup>2</sup>), where h(q) is the one-dimensional Fourier transform of h(x) with q as the wave number, ⟨⟩ denotes the time average, k<sub>B</sub> is the Boltzman constant, and L<sub>x</sub>, L<sub>z</sub> are width and thickness of the simulation box. The interfacial stiffness, γ<sub>m</sub>, is used as a fair estimation of the γ<sub>m</sub>, whereas different orientations of the crystal structure enable the determination of γ<sub>m</sub> anisotropy. Nevertheless, at this point, it is worth mentioning that the direct studies of model<sup>57,59,60,76</sup> and realistic<sup>62,67,77,78</sup> systems suggest that this effect is usually relatively weak, and therefore γ<sub>m</sub> can be obtained from γ̃<sub>m</sub> determined from a single crystallographic orientation. Then, ln(⟨|h(q)|<sup>2</sup>⟩) is a linear function of ln(q) with a slope equal to -2<sup>57</sup>, and the intercept is directly related to γ̃<sub>m</sub>. Thus, γ<sub>m</sub> can be estimated by fitting the obtained dependence of ⟨|h(q)|<sup>2</sup>⟩ on q (expressed in logarithmic scales) to the linear function with the constant slope equal to -2 and analyzing its intercept, see Fig. 2b, where discussed fits are presented for both studied systems. As the initial stage of this analysis we determined the region at which ⟨|h(q)|<sup>2</sup>⟩ is characterized by the linear dependence on ln(q) with slope equal to -2. Subsequently the linear function was adjusted to the chosen points. Estimated in this way values of γ<sub>m</sub> are equal to 2.53 ± 0.44mJ/m<sup>2</sup>
<figure><img src="image_5.png" /><figcaption>**Figure 2.** In the panel (a), the value of the rotational-invariant order parameter along the direction perpendicular to the solid-liquid interface for the system I is presented. In the panel (b), the fluctuation spectrum of the interface height for two RLM systems is shown. The straight lines represent the fit of the linear function with a constant slope equal to -2.</figcaption></figure>

<figure><img src="image_6.png" /><figcaption>**Figure 3.** The values of the nucleation and crystal growth rates for different temperatures are shown. The arrows indicate N and U value obtained at the temperature at which the studied RM system crystallized during cooling.</figcaption></figure>

and 4.15 ± 0.29mJ/m² respectively for the systems I and II. Subsequently, the temperature dependence of γ was estimated according to the Turnbull law⁷⁹,

$$ \gamma(T) = \gamma_m \left( \frac{\rho_{cr}(T)}{\rho_{cr}(T_m)} \right)^{2/3} \left( \frac{\Delta H(T)}{\Delta H_m} \right). $$

Finally, the calculated values of the N and U are presented in Fig. 3.

It can be clearly seen that N and U possess maxima located close to each other for both studied systems. Interestingly, the above maxima are located around the temperature at which the given system crystallizes during cooling. Arrows in Fig. 3 indicate the mentioned temperatures. Thus, CNT fairly predicts the thermodynamic conditions of crystallization. However, at this point, it is worth noting that despite the fact that both systems are very similar and that the procedures of the performed experiments are identical the system II crystallizes at temperatures much lower than $T_m$ ($T_m - T = 45K$ for system I, whereas for system II this difference equals 20K). This observation is even more intriguing when one takes into account that around 20K below the melting temperature, N for system II is about 3 times higher than for system I. Additionally, at discussed thermodynamic conditions U for system II is also much higher than for system I, and therefore from the CNT point of view, neither N nor U suspend the crystallization. Consequently, the system II should easily crystallize much faster than it is observed during cooling experiments. Moreover, we would like to put the reader's attention on another intriguing fact. Examining the cooling procedure, one might observe that despite the fact that at $T_{cr}$ the system
<figure><img src="image_7.png" /><figcaption>**Figure 4.** The time evolution of systems volumes at temperatures 5 K higher (top) and lower (bottom) than the temperature at which RMs systems crystallized during cooling.</figcaption></figure>

II possesses almost 10 times higher maximal values of the N and U, the time needed to crystallization of both systems are similar, i.e., the initial liquid structures become entirely solid within the same simulation time (10 ns).

Discussion

Nevertheless, it must be noted that the crystallization process starts from the stochastic formation of the critical nuclei within the liquid. Therefore, to examine the crystallization tendency in detail, and then to confirm that the characteristics of crystallization process for examining systems are indeed similar we simulate the liquid structure at temperature 5K higher than $T_{cr}$ for 200ns or till the time at which we observed the crystallization event. The results are presented in Fig. 4a, where one can see that increase in the temperature implies an extension of the time needed for registration of the crystallization for both systems.

This observation corresponds with a prediction of the CNT. Additionally, it is worth mentioning that system I is more sensitive for applied temperature changes because 3 from 5 simulation runs did not end in the crystallization, whereas system II always crystallized. Summarizing, we can state that at temperature equals $T_{cr} + 5K$ the time needed for the registration of crystallization is slightly shorter for system II. However, at discussed temperature conditions the nucleation rate is more than 25 times higher for system II than for system I. It implies that assuming that the time needed for crystallization is equal to about 100ns for system II (see Fig. 4a), we could anticipate that system I would persist in liquid state for 2500ns (the total time would be 25 time longer than 100 ns). However, within only 200ns, the system I crystallized twice. Thus, comparing both systems the prediction of CNT does not correspond to the observed results. Subsequently, we simulated both systems at temperature 5K lower than $T_{cr}$. The results are shown in Fig. 4b. As one can observe both systems always crystallize within 5 ns, i.e., within the time which is 2 times shorter than in the case of cooling experiment. Similar to previous results the crystallization process proceeds slightly faster for system II. However, in this case, the differences between both systems are less prominent. Hence, at temperatures 5K lower than $T_{cr}$ the stability behaviors of studied systems can be considered as comparable. Moreover, we would like to note that CNT predicts that at temperature lower than $T_{cr}$ the crystallization process for system II should slow down due to the decrease in N and U, which is not observed in the performed experiments, see Fig. 5b.

At this point it has to be noted that N values presented in Fig. 3 are expressed in the unit of 1/nm³s, which implies that N, and hence the number of the created critical nucleuses depends on the system size. The two studied systems, which are comprised of the same amount of the molecules, are simulated at various temperatures. Therefore, they exhibit different volumes. Nevertheless, as it can be seen in Fig. 4, the differences in v equal only about 5% and therefore its impact on N can be neglected.
<figure><img src="image_8.png" /><figcaption>**Figure 5.** In the left panel the temperature dependences of diffusion constant for RMs systems are presented. The black lines represent the fit of the VFT equation, $log(\frac{1}{D_0}) = log(\frac{1}{D_0}) + log(e)(\frac{B}{T-T_0})$. The fit parameters are $log(\frac{1}{D_0}) = -(1068 \pm 4)10^{-2}$, $B = (52 \pm 3)10^1$K, $T_0 = 35 \pm 4$K for system I and $log(\frac{1}{D_0}) = -(1110 \pm 7)10^{-2}$, $B = (12 \pm 1)10^2$K, $T_0 = 17 \pm 2$K for system II. The right panel shows the values of the nucleation rate divided by the diffusion constant for studied systems.</figcaption></figure>

Putting an attention on Eqs. (1) and (2), the possible explanations of observed differences between prediction of CNT and computational experiment can be found. Both equations consider the diffusion of the molecules. It immediately implies that systems possessing a higher value of *D* exhibit higher values of *N* and *U*. It is crucial in the case of comparison between systems which are characterized by the significant differences in $T_m$ because the diffusion strongly depends on the temperature. The higher $T_m$ implies the faster diffusion of the system and consequently the greater *N* and *U* values are predicted by CNT. The latter seems to be crucial especially in the case of the system characterized by similar structure. In Fig. 5a one can see that $D(T_m)$ for system II is higher for about one decade than for system I.

The latter immediately implies 10 times higher values of *N* and *U* for system II. Interestingly, one can see in Fig. 3 that *N* and *U* between both systems differ also about 10 times. Hence, the reported variation in CNT predictions for both systems is consistent with differences in *D*. Following this observation, in Fig. 5b we present $N/D$ values for both systems, which indeed are very similar. This finding confirms that the main reason for divergences in CNT predictions for studied herein system is the noticeable difference in *D* values. At this point it is worth mentioning that for liquid close to the melting conditions the ratio between rotational and translational diffusion is constant and independent on the temperature⁸⁰. Hence, neither translation nor rotation can be treated as a limiting factor for crystallization process of system I.

Summarizing, in this paper we calculate the *N* and *U* curves according to CNT for two RM systems, which differ exclusively in the value of the dipole moment. The use of proposed model molecules enables entire elimination of the molecular structure role in the crystallization process. Importantly, it makes also that, in contrast to standard simple model systems, obtained results for the system *I* cannot be used to reproduce the results determined for the system *II*. It is also worth noting that, we calculate $\gamma$ using CFM, instead of estimation of its value. Our results show that *N* and *U* curves differ strongly for two studied system. The system with higher value of the dipole moment is characterized by about 10 times higher *N* and *U*. Interestingly, despite the fact that the system *II* exhibits drastically higher values of *N* and *U*, it does not crystallize at expected thermodynamic conditions i.e., at conditions at which *N* and *U* for second system are sufficient to observe the crystallization process. Our results suggest that the main reason for observed discrepancies between results of performed computational experiments and CNT predictions is the diffusion constant.

## Methods

We employ the previously proposed the quasi-real molecules of the rhombus shape, i.e., rhombus-like molecules (RMs), which remarkable advantage is that keeping the simplicity of classical model systems, they display the structural anisotropy typical for the real molecules and simultaneously enable the creation of the differently oriented dipole moments, $\mu^{51,81-83}$. On the basis of the results reported in Ref.⁵¹, we know that only one of the 5 different systems, which vary in the values and orientation of $\mu$, crystallizes. Therefore, we use this system as a reference one. It consists of 4 identical atoms (of carbon atom mass) arranged, as we already mentioned, in a rhombus shape, which implies that RM possess short and long molecular axes (along diagonals of rhombus) simultaneously keeping identical bonds lengths. The latter is set to equal 0.14982 nm, which is close to 0.14 nm, i.e., a length of the bond linking two carbon atoms in the benzene ring. Additionally, the angles between bonds in RM are established to make one diagonal two times longer than the other. To ensure the best mimic of the real molecules by RM, the stiffness of bonds, angles, and dihedrals, as well as the non-bonded interaction between atoms of different RM molecules, are defined by OPLS all-atom force field parameters⁸⁴ provided for carbon atoms of the benzene ring. Then, the permanent $\mu$ is obtained by redefining charges of given atoms, i.e., those arranged along the longer axis are set to 0.0e ($e$ is an elementary charge), whereas those places along shorter one
equal ±0.5e. In this way, we obtain the reference system I (i.e., system C2 from Ref.⁵¹). The second examined herein system, i.e., system II, is identical to the previous one except of the difference in charges' values, which are set to ±0.75e for system II. Consequently, the two model systems differ only in the value of the dipole moment, $\mu_{II} = \frac{3}{2}\mu_I$. Schemes of two studied molecules are presented in the insets of Fig. 1. At this point is worth recalling that consistently to our previous mention the applied method for RM molecules creation makes results obtained for one of them cannot be somehow used to reproduce the results for another one. However, the identity of the molecular structure is preserved.

The heating process was performed by the use of the GROMACS software⁸⁵⁻⁸⁸ at conditions of constant temperature and pressure, which were controlled by the Nose-Hoover thermostat⁸⁹⁻⁹¹ and Martyna-Tucker-man-Tobias-Klein barostat⁹²,⁹³ ($p = 1000$bar) respectively. The increase in the temperature between subsequent steps is equal to 5K. The first half of simulations, which last for 5 ns (the time step $dt = 0.001$ ps), was spared for equilibration of the system, whilst the data was collected during the second half.

Received: 8 December 2021; Accepted: 20 April 2022

Published online: 10 June 2022

## References

1. Yu, L. Amorphous pharmaceutical solids: Preparation, characterization and stabilization. *Adv. Drug Deliv. Rev.* **48**, 27–42 (2001).

2. Demetriou, M. D. et al. A damage-tolerant glass. *Nat. Mater.* **10**, 123–128 (2011).

3. Myerson, A. S. *Handbook of Industrial Crystallization* 2nd edn. (Butterworth Heinemann, 2002).

4. Novoa, J. J., Braga, D. & Addadi, L. *Engineering of Crystalline Materials Properties: State of the Art in Modeling, Design and Applications* (Springer, 2008).

5. Debenedetti, P. G. *Metastable Liquids Concepts and Principles* (Princeton University Press, 1996). https://doi.org/10.1088/0953-8984/15/1/308.

6. Wedekind, J., Strey, R. & Reguera, D. New method to analyze simulations of activated processes. *J. Chem. Phys.* **126**, 134103 (2007).

7. Mokshin, A. V. & Barrat, J.-L. Shear-induced crystallization of an amorphous system. *Phys. Rev. E* **77**, 21505 (2008).

8. Mokshin, A. V. & Barrat, J.-L. Shear induced structural ordering of a model metallic glass. *J. Chem. Phys.* **130**, 34502 (2009).

9. Rein ten Wolde, P., Ruiz-Montero, M. J. & Frenkel, D. Numerical calculation of the rate of crystal nucleation in a Lennard-Jones system at moderate undercooling. *J. Chem. Phys.* **104**, 9932–9947 (1996).

10. Yarullin, D. T., Galimzyanov, B. N. & Mokshin, A. V. Direct evaluation of attachment and detachment rate factors of atoms in crystallizing supercooled liquids. *J. Chem. Phys.* **152**, 224501 (2020).

11. Mokshin, A. V. & Galimzyanov, B. N. Steady-state homogeneous nucleation and growth of water droplets: Extended numerical treatment. *J. Phys. Chem. B* **116**, 11959–11967 (2012).

12. Uhlmann, D. R. Crystallization and melting in glass-forming systems. In *Kinetics of Reactions in Ionic Systems* 172–197 (Springer, 1969). https://doi.org/10.1007/978-1-4899-6461-8_9.

13. Gutzow, I. S. & Schmelzer, J. W. P. States of aggregation, thermodynamic phases, phase transformations, and the vitreous state. *Vitreous State* **20**, 7–67. https://doi.org/10.1007/978-3-642-34633-0_2 (2013).

14. Koperwas, K. et al. Glass-forming tendency of molecular liquids and the strength of the intermolecular attractions. *Sci. Rep.* **6**, 36934 (2016).

15. Weeks, J. D., Chandler, D. & Andersen, H. C. Role of repulsive forces in determining the equilibrium structure of simple liquids. *J. Chem. Phys.* **54**, 5237–5247 (1971).

16. Hsu, C. S. & Rahman, A. Interaction potentials and their effect on crystal nucleation and symmetry. *J. Chem. Phys.* **71**, 4974 (1979).

17. Kob, W. & Andersen, H. C. Testing mode-coupling theory for a supercooled binary Lennard-Jones mixture I: The van Hove correlation function. *Phys. Rev. E* **51**, 4626–4641 (1995).

18. Berthier, L. & Tarjus, G. Nonperturbative effect of attractive forces in viscous liquids. *Phys. Rev. Lett.* **103**, 25 (2009).

19. Baidakov, V. G., Tipeev, A. O., Bobrov, K. S. & Ionov, G. V. Crystal nucleation rate isotherms in Lennard-Jones liquids. *J. Chem. Phys.* **132**, 234505 (2010).

20. Baidakov, V. G. & Tipeev, A. O. Crystal nucleation and the solid-liquid interfacial free energy. *J. Chem. Phys.* **136**, 74510 (2012).

21. Tanaka, K. K., Kawamura, K., Tanaka, H. & Nakazawa, K. Tests of the homogeneous nucleation theory with molecular-dynamics simulations. *I. Lennard-Jones molecules*. *J. Chem. Phys.* **122**, 184514 (2005).

22. Bordat, P., Affouard, F., Descamps, M. & Ngai, K. L. Does the interaction potential determine both the fragility of a liquid and the vibrational properties of its glassy state? *Phys. Rev. Lett.* **93**, 105502 (2004).

23. Shi, Z., Debenedetti, P. G., Stillinger, F. H. & Ginart, P. Structure, dynamics, and thermodynamics of a family of potentials with tunable softness. *J. Chem. Phys.* **135**, 084513 (2011).

24. Zhang, K. et al. Computational studies of the glass-forming ability of model bulk metallic glasses. *J. Chem. Phys.* **139**, 124503 (2013).

25. Koperwas, K. et al. Influence of pressure on the crystallization of systems characterized by different intermolecular attraction. *Phys. Rev. B* **96**, 224106 (2017).

26. Toxvaerd, S. Role of the attractive forces in a supercooled liquid. *Phys. Rev. E* **103**, 22611 (2021).

27. Toxvaerd, S. The role of local bond-order at crystallization in a simple supercooled liquid. *Eur. Phys. J. B* **93**, 202 (2020).

28. Singh, A. & Singh, Y. How attractive and repulsive interactions affect structure ordering and dynamics of glass-forming liquids. *Phys. Rev. E* **103**, 52105 (2021).

29. Auer, S. & Frenkel, D. Crystallization of weakly charged colloidal spheres: A numerical study. *J. Phys. Condens. Matter* **14**, 7667 (2002).

30. Valdès, L.-C., Gerges, J., Mizuguchi, T. & Affouard, F. Crystallization tendencies of modelled Lennard-Jones liquids with different attractions. *J. Chem. Phys.* **148**, 014501 (2018).

31. Vega, C., Paras, E. P. A. & Monson, P. A. Solid-fluid equilibria for hard dumbbells via Monte Carlo simulation. *J. Chem. Phys.* **96**, 9060–9072 (1992).

32. Vega, C., Paras, E. P. A. & Monson, P. A. On the stability of the plastic crystal phase of hard dumbbell solids. *J. Chem. Phys.* **97**, 8543–8548 (1992).

33. Vega, C. & Monson, P. A. Plastic crystal phases of hard dumbbells and hard spherocylinders. *J. Chem. Phys.* **107**, 2696–2697 (1997).

34. Marechal, M. & Dijkstra, M. Stability of orientationally disordered crystal structures of colloidal hard dumbbells. *Phys. Rev. E* **77**, 061405 (2008).

35. Zubieta, I., Vázquez del Saz, M., Llombart, P., Vega, C. & Noya, E. G. Nucleation of pseudo hard-spheres and dumbbells at moderate metastability: Appearance of A15 Frank-Kasper phase at intermediate elongations. *Phys. Chem. Chem. Phys.* **21**, 1656–1670 (2019).

36. Bolhuis, P. & Frenkel, D. Tracing the phase boundaries of hard spherocylinders. *J. Chem. Phys.* **106**, 666–687 (1997).
37. de Miguel, E., del Rio, E., Brown, J. T. & Allen, M. P. Effect of the attractive interactions on the phase behavior of the Gay-Berne liquid crystal model. *J. Chem. Phys.* **105**, 4234–4249 (1996). https://doi.org/10.1002/9780470143513.ch3.

38. Kihara, T. Convex molecules in gaseous and crystalline states. In *Advances in Chemical Physics* (ed. Prigogine, I.) 147–188 (Wiley, 2007). https://doi.org/10.1002/9780470143513.ch3.

39. Berne, B. J. & Pechukas, P. Gaussian model potentials for molecular interactions. *J. Chem. Phys.* **56**, 4213–4216 (1972).

40. Gay, J. G. & Berne, B. J. Modification of the overlap potential to mimic a linear site–site potential. *J. Chem. Phys.* **74**, 3316–3319 (1981).

41. Thompson, S. M. & Gubbins, K. E. Structure of the liquid–vapor interface of molecular fluids: Homonuclear diatomic molecules. *J. Chem. Phys.* **74**, 6467–6479 (1981).

42. Thompson, S. M. & Gubbins, K. E. Molecular orientation at a vapor–liquid interface: Theoretical and computer simulation results for a model of chlorine. *J. Chem. Phys.* **70**, 4947–4951 (1979).

43. Nijmeijer, M. J. P., Bakker, A. F., Bruin, C. & Sikkenk, J. H. A molecular dynamics simulation of the Lennard-Jones liquid–vapor interface. *J. Chem. Phys.* **89**, 3789–3792 (1988).

44. Rao, M. & Berne, B. J. On the location of surface of tension in the planar interface between liquid and vapour. *Mol. Phys.* **37**, 455–461 (1979).

45. Cheung, P. S. Y. & Powles, J. G. The properties of liquid nitrogen. *Mol. Phys.* **30**, 921–949 (1975).

46. Singer, K., Taylor, A. & Singer, J. V. L. Thermodynamic and structural properties of liquids modelled by ‘2-Lennard-Jones centres’ pair potentials. *Mol. Phys.* **33**, 1757–1795 (1977).

47. Romano, S. & Singer, K. Calculation of the entropy of liquid chlorine and bromine by computer simulation. *Mol. Phys.* **37**, 1765–1772 (1979).

48. Monson, P. A. Thermodynamic properties of molecular fluids from the site-site Ornstein-Zernike equation. *Mol. Phys.* **53**, 1209–1223 (1984).

49. McGuigan, D. B., Lupkowski, M., Paquet, D. M. & Monson, P. A. Phase diagrams of interaction site fluids. *Mol. Phys.* **67**, 33–52 (1989).

50. Rull, L. F. Phase diagram of a liquid crystal model: A computer simulation study. *Phys. A Stat. Mech. Appl.* **220**, 113–138 (1995).

51. Koperwas, K., Adrjanowicz, K., Grzybowski, A. & Paluch, M. The role of the dipole moment orientations in the crystallization tendency of the van der Waals liquids—molecular dynamics simulations. *Sci. Rep.* **10**, 283 (2020).

52. Gutzow, I., Durschang, B. & Rüssel, C. Crystallization of glassforming melts under hydrostatic pressure and shear stress: Part I Crystallization catalysis under hydrostatic pressure: Possibilities and limitations. *J. Mater. Sci.* **32**, 5389–5403 (1997).

53. Adrjanowicz, K., Grzybowski, A., Grzybowska, K., Pionteck, J. & Paluch, M. Effect of High pressure on crystallization kinetics of Van der Waals liquid: An experimental and theoretical study. *Cryst. Growth Des.* **14**, 2097–2104 (2014).

54. Broughton, J. Q. & Gilmer, G. H. Molecular dynamics investigation of the crystal–fluid interface. III. Dynamical properties of fcc crystal–vapor systems. *J. Chem. Phys.* **79**, 5119–5127 (1983).

55. Davidchack, R. L. & Laird, B. B. Direct calculation of the hard-sphere crystal melt interfacial free energy. *Phys. Rev. Lett.* **85**, 4751–4754 (2000).

56. Davidchack, R. L. & Laird, B. B. Direct calculation of the crystal–melt interfacial free energies for continuous potentials: Application to the Lennard-Jones system. *J. Chem. Phys.* **118**, 7651 (2003).

57. Hoyt, J. J., Asta, M. & Karma, A. Method for computing the anisotropy of the solid-liquid interfacial free energy. *Phys. Rev. Lett.* **86**, 5530–5533 (2001).

58. Hoyt, J. J. & Asta, M. Atomistic computation of liquid diffusivity, solid-liquid interfacial free energy, and kinetic coefficient in Au and Ag. *Phys. Rev. B* **65**, 214106 (2002).

59. Morris, J. R. Complete mapping of the anisotropic free energy of the crystal-melt interface in Al. *Phys. Rev. B* **66**, 144104 (2002).

60. Morris, J. R. & Song, X. The anisotropic free energy of the Lennard-Jones crystal-melt interface. *J. Chem. Phys.* **119**, 3920–3925 (2003).

61. Mu, Y., Houk, A. & Song, X. Anisotropic interfacial free energies of the hard-sphere crystal–melt interfaces†. *J. Phys. Chem. B* **109**, 6500–6504 (2005).

62. Handel, R., Davidchack, R. L., Anwar, J. & Brukhno, A. Direct calculation of solid-liquid interfacial free energy for molecular systems: TIP4P ice-water interface. *Phys. Rev. Lett.* **100**, 036104 (2008).

63. Ramsteiner, I. B., Weitz, D. A. & Spaepen, F. Stiffness of the crystal-liquid interface in a hard-sphere colloidal system measured from capillary fluctuations. *Phys. Rev. E* **82**, 41603 (2010).

64. Rozas, R. E. & Horbach, J. Capillary wave analysis of rough solid-liquid interfaces in nickel. *Europhys. Lett.* **93**, 26006 (2011).

65. Asta, M., Hoyt, J. J. & Karma, A. Calculation of alloy solid-liquid interfacial free energies from atomic-scale simulations. *Phys. Rev. B* **66**, 100101 (2002).

66. Becker, C. A., Olmsted, D. L., Asta, M., Hoyt, J. J. & Foiles, S. M. Atomistic simulations of crystal-melt interfaces in a model binary alloy: Interfacial free energies, adsorption coefficients, and excess entropy. *Phys. Rev. B* **79**, 54109 (2009).

67. Feng, X. & Laird, B. B. Calculation of the crystal-melt interfacial free energy of succinonitrile from molecular simulation. *J. Chem. Phys.* **124**, 044707 (2006).

68. Gerges, J. & Affouard, F. Predictive calculation of the crystallization tendency of model pharmaceuticals in the supercooled state from molecular dynamics simulations. *J. Phys. Chem. B* **119**, 10768–10783 (2015).

69. Gerges, J. & Affouard, F. Insight from molecular dynamics simulations on the crystallization tendency of indomethacin polymorphs in the undercooled liquid state. *J. Pharm. Sci.* **109**, 1086–1095 (2020).

70. Koperwas, K. et al. Pressure dependence of the crystallization rate for the S-enantiomer and a racemic mixture of ibuprofen. *Cryst. Growth Des.* https://doi.org/10.1021/acs.cgd.1c00980 (2021).

71. Steinhardt, P. J., Nelson, D. R. & Ronchetti, M. Bond-orientational order in liquids and glasses. *Phys. Rev. B* **28**, 784–805 (1983).

72. Chushak, Y. & Bartell, L. S. Crystal nucleation and growth in large clusters of SeF6 from molecular dynamics simulations. *J. Phys. Chem. A* **104**, 9328–9336 (2000).

73. Lechner, W. & Dellago, C. Accurate determination of crystal structures based on averaged local bond order parameters. *J. Chem. Phys.* **129**, 114707 (2008).

74. Reinhardt, A., Doye, J. P. K., Noya, E. G. & Vega, C. Local order parameters for use in driving homogeneous ice nucleation with all-atom models of water. *J. Chem. Phys.* **137**, 194504 (2012).

75. Mickel, W., Kapfer, S. C., Schröder-Turk, G. E. & Mecke, K. Shortcomings of the bond orientational order parameters for the analysis of disordered particulate matter. *J. Chem. Phys.* **138**, 44501 (2013).

76. Amini, M. & Laird, B. B. Crystal-melt interfacial free energy of binary hard spheres from capillary fluctuations. *Phys. Rev. B* **78**, 144112 (2008).

77. Benet, J., MacDowell, L. G. & Sanz, E. A study of the ice-water interface using the TIP4P/2005 water model. *Phys. Chem. Chem. Phys.* **16**, 22159–22166 (2014).

78. Davidchack, R. L., Handel, R., Anwar, J. & Brukhno, A. V. Ice I h–water interfacial free energy of simple water Models with Full Electrostatic Interactions. *J. Chem. Theory Comput.* **8**, 2383–2390 (2012).

79. Turnbull, D. Formation of crystal nuclei in liquid metals. *J. Appl. Phys.* **21**, 1022–1028 (1950).

80. Edmond, K. V., Elsesser, M. T., Hunter, G. L., Pine, D. J. & Weeks, E. R. Decoupling of rotational and translational diffusion in supercooled colloidal fluids. *Proc. Natl. Acad. Sci.* **109**, 17891–17896 (2012).
81. Koperwas, K., Grzybowski, A. & Paluch, M. The effect of molecular architecture on the physical properties of supercooled liquids studied by MD simulations: Density scaling and its relation to the equation of state. *J. Chem. Phys.* **150**, 014501 (2019).

82. Koperwas, K., Grzybowski, A. & Paluch, M. Exploring the connection between the density-scaling exponent and the intermolecular potential for liquids on the basis of computer simulations of quasireal model systems. *Phys. Rev. E* **101**, 012613 (2020).

83. Koperwas, K., Grzybowski, A. & Paluch, M. Virial-potential-energy correlation and its relation to density scaling for quasireal model systems. *Phys. Rev. E* **102**, 62140 (2020).

84. Jorgensen, W. L., Maxwell, D. S. & Tirado-Rives, J. Development and testing of the OPLS all-atom force field on conformational energetics and properties of organic liquids. *J. Am. Chem. Soc.* **118**, 11225-11236 (1996).

85. Berendsen, H. J. C., van der Spoel, D. & van Drunen, R. GROMACS: A message-passing parallel molecular dynamics implementation. *Comput. Phys. Commun.* **91**, 43–56 (1995).

86. Van Der Spoel, D. et al. GROMACS: Fast, flexible, and free. *J. Comput. Chem.* https://doi.org/10.1002/jcc.20291 (2005).

87. Hess, B., Kutzner, C. van der Spoel, D. & Lindahl, E. GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation. *J. Chem. Theory Comput.* **4**, 435–447 (2008).

88. Pronk, S. et al. GROMACS 4.5: A high-throughput and highly parallel open source molecular simulation toolkit. *Bioinformatics* **29**, 845–854 (2013).

89. Nosé, S. A molecular dynamics method for simulations in the canonical ensemble. *Mol. Phys.* **52**, 255–268 (1984).

90. Nosé, S. A unified formulation of the constant temperature molecular dynamics methods. *J. Chem. Phys.* **81**, 511–519 (1984).

91. Hoover, W. G. Canonical dynamics: Equilibrium phase-space distributions. *Phys. Rev. A* **31**, 1695–1697 (1985).

92. Martyna, G. J., Tobias, D. J. & Klein, M. L. Constant pressure molecular dynamics algorithms. *J. Chem. Phys.* **101**, 4177–4189 (1994).

93. Martyna, G. J., Tuckerman, M. E., Tobias, D. J. & Klein, M. L. Explicit reversible integrators for extended systems dynamics. *Mol. Phys.* **87**, 1117–1157 (1996).

## Acknowledgements

K.K., F.K., A.G., and M.P. are deeply grateful for the financial support by the Polish National Science Centre within the framework of the Maestro10 project (Grant No. UMO-2018/30/A/ST3/00323). This project has also received funding from the Interreg 2 Seas program 2014-2020 co-funded by the European Regional Development Fund (FEDER) under subsidy contract 2S01-059_IMODE.

## Author contributions

K.K. developed the concept of research and drafted the manuscript. F.K. and K.K. carried out the computations, which were supervised by F.A., whereas A.G. and M.P. directed entire project. All authors read and approved the final manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

Correspondence and requests for materials should be addressed to K.K. or F.K.

Reprints and permissions information is available at www.nature.com/reprints.

Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_9.png" /></figure>

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2022