Hard spheres and hard disks are widely accepted as the first choice approximation to model a variety of soft condensed matter objects [1]. In spite of being a rather simple representation of the substance, hard core-based model systems still are capable to recover a number of basic properties and effects related to the structure, thermodynamics (e.g., phase transitions), and dynamics of real systems. Recently, the interest has been revived in the properties of hard sphere fluid confined to a narrow channel of the width that does not exceed two hard core diameters. In such system, commonly referred to as quasi-one dimensional (quasi-1D) system [2], the hard-core particles cannot pass the nearest neighbors and their motion is restricted by the neighbors. Of a particular interest is the so-called single-file quasi-1D system with the width that does not exceed (1+3/2) of disk diameter [3] as then a disk cannot touch more than one neighbor from each side. This is a substantial simplification that allows one to make a contact with the exact Tonks solution for the purely 1D hard rod system [4]. The reasons for the interest to this system are both basic and applied. The fundamental interest stems from the existence of the analytical transfer matrix approach for the isobaric partition function [5–9] and the exact canonical partition function of this system, which is recently derived in [10]. At present, the theoretical research in this area has been mostly concerned with the efforts (i) to get an insight into the mechanism that governs transformation of the properties of hard core systems as their dimensionality approaches 2 and 3 [5, 11], and (ii) to consider these systems as glass formers, see [7–9, 12–14] and review [15]. The practical interest stems from the possibility to use such a simple model to capture properties of more complex systems, e.g., to explain diffusion in zeolite and carbon channels [16–18], microfluidic devices [19], in the technology of bio-integrated nanodevices [3] etc. by treating the finite length axis of the quasi-1D system as the pore width.

Confined many-particle systems are very versatile and complex, and computer simulations are a perfect and in the most cases only available tool to study their structure and dynamics. It has been widely applied to study the structure and single-particle dynamics in quasi-1D systems with various model smooth/continuous particle–particle and particle–wall interaction potentials. But even the quasi-1D systems partially involving hard core particles and the quasi-1D systems of pure hard core particles (that are the main subject of present study) can be so different in their physics that comparison of apparently similar effects in such systems often has no much physical sense. For instance, in a quasi-1D system of hard core particles suspended in a viscoelastic liquid solvent their effective interactions are mainly of a hydrodynamic nature and the role of confining walls is not so much in restricting particles' motion as in setting the boundary condition on the liquid solvent's flow. At the same time, in a quasi-1D system of pure hard core particles all interactions are of an entropic nature. Another example concerns a zigzag particle arrangement that very often is the object of fundamental interest for different quasi-1D systems. However, because of the difference in interaction potentials, the general idea of the zigzag geometry is often the only thing in common. For instance, a zigzag structure can occur both in quasi-1D system with extremely short-range interaction such as a hard core repulsion (e.g., see [9]) and in quasi-1D system with long-range interaction such as a screened electrostatic repulsion in dusty plasma (e.g., see review [20]). The difference between these two interactions is known to be fundamental even in a purely 1D geometry in which the phase transition exists only for a long-range particle–particle repulsion [21]. At the same time, the computer simulations studies of dynamic properties of the quasi-1D systems that contain only hard core particles, i.e., the quasi-1D hard core systems, have been primarily performed in the context of glassy dynamics. They were mostly focused on the problems specific for the single particle dynamics such as disk hopping time between different quasi-equilibrium states and defects' dynamics [7, 9, 12], and time dependence of the particle displacement [9, 13]. The glassy behavior of hard disk fluid confined to a narrow hard wall channel was also studied theoretically in [7–9] by analyzing the transverse and longitudinal equilibrium static pair correlation functions by the transfer matrix approach [5].

Our recent computer simulations studies [22, 23] have addressed the collective time correlation functions in the bulk 2D and 3D hard core systems. What in the following is referred to as collective dynamics means correlations in the cooperative motion of a many-particle system. Such an approach to understanding the dynamical processes is general for studies of propagating waves (sound, shear, and heat ones) and relaxation processes (thermal relaxation, structural relaxation, stress relaxation, etc.). The corresponding collective dynamic variables are defined via fluctuations of the conserved quantities: number of particles, longitudinal and transverse components of total momentum, and energy. All these collective dynamic variables are known from the hydrodynamic approach and can be used for theoretical description of the long-wavelength processes. For the theoretical description of collective dynamics beyond the hydrodynamic regime, the set of collective variables is extended by the orthogonal ones (orthogonal to hydrodynamic variables), which represent the longitudinal and transverse components of stress tensor, energy current, etc. More precisely, computer simulations [22, 23] were used to find how collective dynamics of bulk hard spheres and bulk hard disks behaves on different spatial scales. Rather unexpectedly the short-wavelength shear waves were found in both cases from the well-defined peaks of the transverse current spectral functions. Although nature of shear waves in hard-core systems appears to be essentially the same as in simple fluids [24], i.e., they emerge as the short-wavelength excitations due to coupling of the transverse current and transverse component of stress tensor, the hard-core fluids are known for the absence of viscoelastic effects. Interestingly, while there are no short-wavelength shear waves in a 2D hard disk fluid at low particle densities, these were observed in the range of higher densities by showing certain particular features just before the freezing transition.

The present paper is devoted to a quasi-1D hard disk system. An essential difference between the methods to study collective dynamics in hard-core systems and those for models based on analytic/continuous two-body potentials stems from the absence of a local energy minima in the former case. Even in the case of solids composed of hard spheres or hard disks, the particles do not oscillate around the minimum of potential well, but are moving ballistically in the cage formed by the nearest neighbors. This difference gives rise to an interesting specific aspect of the collective dynamics in a hard core system, namely, an existence of strong correlations between the emergence of short-wavelength shear waves and the caging phenomenon [25, 26], which was noticed in the case of 2D hard disk system [22]. As caging inevitably emerges in a hard disk system under quasi-1D confinement, it is quite natural to expect the existence of short-wavelength shear waves in this case as well. Then it is not clear how the collective modes behave when additional confinement is imposed, how the reflections from channel boundaries affect the longitudinal and transverse excitations originated from particle–particle collisions, how the single-particle ballistic motions sum up to form a collective oscillation modes. The case of transverse excitations is of particular interest since latter are not present in the 1D prototype of the quasi-1D hard disk system. All these issues remain unexplored. Therefore, our aim is to perform molecular dynamics simulations of a quasi-1D system of hard disks, calculate its static structural properties, and make a link with the collective dynamics of the system. The rest of the paper is organized as follows. In the next section, we present the information on molecular dynamics simulations; section 3 contains results of the static and dynamic properties of the studied systems. In the last section, we discuss our findings.

The quasi-1D system is modeled by placing N hard disks of diameter σ into elongated rectangular box formed by two walls (lines) of length L_x ≡ L that are separated by a distance L_y ≡ H = σ + h with h < σ such that disks cannot pass each other. The disk-disk two-body interaction is given by

where r_ij = |r_j − r_i| is the distance between disk centers. Additionally, both confining walls are impenetrable,

thereby forming the hard-wall channel (or pore) of the width H = σ + h and of the length L. The ends of the channel are open and the periodic boundary conditions in x-direction are employed.

Taking into account that width H = σ corresponds to 1D case, the range of channel widths σ < H < 2σ could be considered as a bridge between 1D and higher dimensions. Among continuous variety of the quasi-1D system width in this range σ < H < 2σ there are two, H/σ = 1.5 and H/σ=1+3/2, when disk ordering at close packing is commensurate with 2D triangular lattice like milestones on the way from 1D to 2D. To illustrate, Figure 1 presents two triangular hard disk arrays that correspond to the most common horizontal (Figure 1A) and vertical (Figure 1B) orientations of 2D triangular lattice that differ by an angle of 30 degrees. The quasi-1D systems that correspond to H/σ = 1.5 and H/σ=1+3/2 are shown by hollow disks that at close packing form two very distinct crystalline zigzags. The main distinction concerns a number of nearest neighbors interacting with each disk, i.e., two for H/σ = 1.5 and four in the case of H/σ=1+3/2. The narrower quasi-1D systems is more close to the 1D system, while wider quasi-1D system is more close to the 2D system. In the present study, we are considering the narrower quasi-1D system with width fixed at H/σ = 1.5.

The main body of simulations, in particular, those that concern the collective dynamics, were performed by using a collection of N = 200 hard disks at the fixed channel width H with h/σ = 0.5 and the varied length, L/σ = 400, 350, 300, 250, 220, 198, 190, and 180. Thus, eight quasi-1D hard disk systems that are characterized by different disk linear density,

were simulated. For convenience of comparison and discussion, particular values of the linear density l, that corresponds to each channel length L, are shown in Table 1 together with corresponding values of the disk number density ρ = Nσ²/(HL) and packing fraction η = Nπσ²/(4HL), since the latter are often used in the literature on quasi-1D systems. The simulation runs with larger systems up to N = 2, 000 were performed as well and are properly indicated in the text. Throughout the paper, the disk hard core diameter σ is used as unit of length, while time is in units of (βmσ²)^1/2. In simulations we used β = σ = m = 1, where β = 1/kT and m is the mass of a disk.

Initially chosen N disks were located randomly inside the channel of the shortest length L, i.e., the highest density considered. Then initial disk configurations for lower densities were obtained by increasing the channel length L. To handle the collisions of hard disks with each other and with the channel hard walls, we employ the event driven molecular dynamics (MD) algorithm [27, 28]. According to this technique, the temperature is kept constant by scaling appropriately the magnitude of velocities of each hard disk such that kinetic energy of the system agrees with the equipartition theorem. The directions of the velocities of disk particles at the beginning of each run were chosen randomly. Before calculating the averages for static quantities, we ran MD simulations with assigned velocities, as described before, to equilibrate the system.

Collective dynamics was studied via calculation of the density-density, energy-energy, and longitudinal (L) and transverse (T) current–current time correlation functions. In order to get equally time-sampled positions of hard disks, i.e., the disk trajectories, the positions of hard disks between the collisional events were interpolated. Proceeding in this way, for the systems of N = 200, 400, and 1,000 hard disks, we recorded the trajectories and velocities along trajectory for each hard disk, totally having dumped 100, 000 configurations for the systems of N = 200 particles, but 40,000 configurations for the systems of N = 400 particles, and only 20,000 configurations for the systems of N = 1, 000 particles with the fixed reduced time interval of 0.01.

The following Fourier components were sampled for each configuration. The Fourier component of the particle density

of the longitudinal component of current density

and the transverse component of current density

as well as of the energy density

Here x_i is the position of disk i along the channel, v_{x, i}(t), v_{y, i}(t) are the components of the ith disk velocity along the channel and perpendicular to the channel, respectively, while εikin is the kinetic energy of ith disk at time t. The wave vector k is defined along the channel x−axis, compatible with the periodicity of simulation box, as k ≡ k_x = 2πm/L_x with m = 1, 2, 3, …. Note that, due to a non-zero channel width H, the positions and velocities of particles in quasi-1D systems have two components, along (x) and perpendicular (y) to the channel. However, we are using only wave vectors sampled along the channel. The reason is that transverse collective excitations in atomistic systems can propagate if at least two particles take part in the transverse collective motion, i.e., the wave number is below π/σ. In the case of a narrow channel of width H = 1.5σ the smallest wavenumber, k_y = 2π/H, is too large.

Having the Fourier components of energy and particle densities, it is straightforward to calculate the Fourier components of the heat density [24]

Here f_ne(k) and f_nn(k) ≡ S(k) are the static energy–density and density–density correlators, respectively, while S(k) is the static structure factor. The fluctuations of heat density permit calculations of wavenumber-dependent specific heat at constant volume C_v(k), which in the long-wavelength limit tends to its macroscopic value. Another important quantity that can be obtained from observed heat density dynamics is the Landau–Placzek–like ratio, which gives information on the share of contributions from relaxing and propagating processes to specific heat C_v [24, 29].

The above-defined time-dependent Fourier components of corresponding densities, equations 4 - 8, describe fluctuations of conserved quantities in monoatomic fluids and form the set of hydrodynamic variables of macroscopic collective dynamics. Having dynamic variables (4–8), we calculated the time-dependent density–density correlation function,

time dependent longitudinal and transverse current–current correlation functions,

and time dependent heat density autocorrelation function,

In what follows, the correlation functions (9–11) were used for numerical time-Fourier transformation in order to obtain the density–density dynamic structure factors S(k, ω), longitudinal/transverse current spectral functions, C^L/T(k, ω) and the heat-density dynamic structure factor S_hh(k, ω).

The system of hard disks in a narrow hard wall channel is anisotropic [5]. The anisotropy is expected due to the system setup L >> H, which results in two different pressure components, i.e., the longitudinal, P_L = F_L/H, and transverse, P_T = F_T/L. Here, F_L and F_T are the force per unit cross-section exerted along the channel length L, and the force on a segment of the horizontal wall of length L/(Nσ) = 1/l, respectively. These forces are both of entropic origin and rather sensitive to be evaluated from computer simulations. Fortunately, these forces can be found from the analytical canonical partition function of a quasi-1D hard disk system reported in [10].

Figure 2 presents the dependencies P_L(l) and P_T(l), as well as F_L(l) and F_T(l) for three channel widths, H = σ + h, with h fixed at: (a) h/σ = 0.141, close to the 1D case; (b) h/σ = 0.5, which is far from both lower 1D and higher single-file limits; and (c) h/σ = 0.866, which is very close to the terminate width h/σ=3/2~0.866025 of a single-file system. At low linear density, the transverse force F_T is density independent, implying that the vertical disks' motion is ballistic, the vertical free path is maximum, i.e., H − σ, the disks bounce between channel walls and do not collide with other disks. As for the same low linear density l ≲ 0.8, the transverse force is lower for a wider channel, the frequency of disk bouncing between the horizontal walls has to be lower for wider channels as well.

We see that, for low densities, the transverse pressure P_T is higher than the one along the channel. This is because in that case P_L is determined by a large disk separation along the channel, whereas P_T is determined by a short range of transverse motion bounded from above by H − σ. As linear density increases toward l ≥ 1, disks start to mount one upon another, the vertical free paths decrease and transverse pressure and force both rapidly increase. For h/σ = 0.866, Figure 2C, this results in P_T to be always higher than P_L. However, for more narrow channels with h/σ = 0.141 and h/σ = 0.5, Figures 2A,B, before this happens, at certain value of linear density, which is different for each channel width H, the gaps between disks and walls in transverse direction and the gaps between neighbor disks along the channel become equal, and the pressures on the vertical and horizontal boundaries coincide.

The density distribution profiles, n_y(y), in the transverse y−direction, obtained both from the analytical partition function [10] and from MD simulations for quasi-1D hard disk system of the width H/σ = 1.5, are shown in Figure 3. While linear density is low, l ≲ 0.8, the system is roughly homogeneous across the channel as the density profiles n_y(y) are nearly constant and equal to the correspondent linear density l. In contrast, as linear density increases, l > 1, the disk distribution across the channel shows the tendency of increase in the regions close to the channel walls and decrease in the middle region of the channel. For the highest studied linear density, l = 1.111, the density profile n_y(y) exhibits an almost δ−like shape in the close proximity of the channel walls and practically vanishes elsewhere.

The longitudinal static structure factor, S(k), is calculated from MD simulations in a standard way as instantaneous-time density–density correlator,

In Figure 4, we show the changes in the first peak of S(k) with decrease of linear density l. Regarding the collective dynamics, one of the most important features is location K_max of the main peak of S(k), since the value of k = K_max has the meaning of a pseudo-Brillouin zone boundary in the considered quasi-1D hard disk fluid at a particular density. In this region of wavenumbers, de Gennes's slowing down of the density fluctuations takes place, which is ultimately reflected in the long tails of the density–density time correlation functions. One can see that for linear densities l > 1 the structure factor S(k) is typical for distorted crystals, with the main peak shaped as the sheared-out delta function. For linear densities l < 1, one observes typical fluid-like structure factor. In disordered systems where a structural transition takes place, one can find different slopes of the main peak of S(k) on both sides of the transition [30]. In Figure 5, we show the main peak position K_max of longitudinal static structure factor S(k) as a function of the linear size L of system that is proportional to the inverse 1/l of linear density. The location of the main peak is changing from K_max ≈ 7 for the case of L/σ = 180 (the highest linear density l = 1.111) down to the value K_max ≈ 4.5 for L/σ = 400 (the lowest considered linear density l = 0.5). Indeed, one can see that there is a kink in the behavior of the maxima positions, which occurs in the region of linear densities 0.8 < l < 1.01. According to Figure 2B, this range of linear densities corresponds in our system to the thermodynamic state, where the transverse pressure P_T is lower than the longitudinal pressure P_L. It worth to note that in a quasi-1D system with channel width H/σ = 1.866 the latter newer happens, i.e., always P_T > P_L (see Figure 2C), while there is range of linear densities, 0.93 ≲ l ≲ 1.01, where transverse force F_T is slightly smaller than longitudinal force F_L. This density range corresponds to packing fraction range from η = 0.35 to 0.45 (see Figure 10 of [9]), where the time dependence of the mean-square displacement starts to develop a power-law time dependence due to the slow diffusion of defects in the zigzag arrangement of disks, suggested in [12] as a mechanism for the α relaxation.

The single-particle dynamics can be studied by means of the velocity autocorrelation functions

The velocity autocorrelation functions were already obtained by means of molecular dynamics simulations for the bulk 3D hard sphere [31] and 2D hard disk [22] fluids. In both cases, the authors found that there exists the packing fraction value (η = 0.45 for 3D case and η = 0.65 for 2D hard disks) above which the velocity autocorrelation function at short times develops a negative minimum signaling of a nascent caging. In our case, an anisotropic quasi-1D hard disk system shows essentially different behavior for the directions along and across the channel. In Figure 6, we show the velocity autocorrelation function ψ(t), its xx- and yy-components, and how they depend on the disk density. The xx-component, ψ_xx(t), in the quasi-1D hard disk system behaves in a way very similar to that for bulk 2D hard disk system [22]. It has long-time tails for small linear densities l ≲ 0.667 and changes to a shallow negative minimum at the largest studied linear density l = 1.111 (packing fraction η = 0.582) due to collisions with the nearest neighbors.

The effect of reflections from the channel hard walls is well seen in the transverse component, ψ_yy(t), of the velocity autocorrelation function. It took us by surprise that the characteristic oscillation of ψ_yy(t) due to the wall reflections changed its shape, and especially so for the intermediate density l = 1.01. Toward shorter times, it became more shallow, transforming for the most dense system, l = 1.111, into a very high-frequency heavily damped oscillation. This effect is much better seen in the Fourier-spectrum Z~yy(ω) of the function ψ_yy(t) shown in Figure 7. The characteristic oscillation due to hard wall reflection is changed starting from l = 1.01 to shift toward higher frequencies with increasing oscillation damping. This means that at high densities the zigzag structuring of hard disks prevents their reflections from both channel walls. Instead under zigzag ordering hard disks are reflected from the single nearest wall and from two nearest zigzag neighbors.

The xx- and yy-components of the velocity autocorrelation function in a quasi-1D system were discussed so far only for the case of hard core particles suspended in a viscoelastic liquid solvent [32]. Similar to our Figure 6B, the authors observed a negative minimum and a negative long time tail of the asymptotic form ~ −t^−3/2 for the component parallel to the channel axis. They also concluded that both findings take place only for sticky boundary condition and does not occur if the solvent can slip over the walls. Our finding of a negative minimum of the velocity autocorrelation function in Figure 6B at the highest studied linear density l = 1.111 is of a completely different nature as our system consists solely of hard disks, which fly ballistically between hitting hard obstacles. We did try to estimate the exponentials for a negative long-time tail of ψ_xx(t). However, even in our case of saved 100,000 configurations with N = 200 particles, the noise in the tails was very strong, which did not allow us to reliably estimate the exponents. This shows that the lattice Boltzmann simulations [32] are more appropriate tool for this purpose than the Newtonian MD simulations.

Collective dynamics is usually studied via analysis of the time correlation functions, which store the entire information about collective excitations in the system and their coupling. So far, the collective dynamics has been well-understood on the macroscopic scales in bulk systems, but practically no information is available in the literature on collective excitations in confined low-dimensional systems consisting purely of hard-core particles. The time correlations in the latter case can be essentially different from that in the bulk systems. A simple example of this difference is seen in the transverse dynamics of considered quasi-1D hard disks system, which at the smallest wavenumber accessible in our simulations always shows damped oscillations due to reflection from the channel hard walls. In contrast, in bulk 3d hard-sphere and 2D hard-disk systems, where macroscopic hydrodynamics is valid, no transverse excitations exist at small wavenumbers accessible in simulations. It is seen in Figure 8 that the change of collective transverse current autocorrelation functions with disk density is practically the same as for the single-particle yy-velocity correlations. This makes it evidence that the leading contribution to the transverse current–current time correlation functions comes from the particle reflections from the channel hard walls.

The longitudinal current time correlation functions allow one to estimate the speed of sound for the smallest wavenumbers accessible in simulations. In Figure 9, we show the obtained dependence of the speed of sound on the channel length L that is proportional to inverse 1/l of linear density (see Table 1). As expected, the speed of sound shows a monotonic decrease with increasing channel length, i.e., with decreasing density.

Dispersion of the longitudinal and transverse excitations in the studied quasi-1D hard disk system were obtained from the peak positions of the longitudinal and transverse current spectral functions C^L/T(k, ω) that are the time-Fourier transforms of the MD-derived longitudinal/transverse current–current time correlation functions FJJL/T(k,t). By employing well-established methodology [24], the latter were analyzed for their peak locations, which for different wave numbers k define the dispersion, ω_L/T(k), of longitudinal and transverse excitations. Typical shapes of the spectral function C^L/T(k, ω) are shown in Figure 10 for the system of N = 200 particles and linear density l = 1.01. As the shapes of the spectral functions are noisy, to locate their peak positions and maxima, we made use of the standard Bezier fit for noisy data. The simulations with larger numbers of particles allowed us to access smaller wavenumbers while the dispersion relations within the error bars remained the same.

Figure 11 show the dispersions ω_L/T(k) of longitudinal and transverse excitations in the quasi-1D hard disk system. The parts of Figures 11A,B present data for the linear densities l < 1. At first glance in this case the dispersion ω_L(k) of longitudinal excitations in quasi-1D hard disk system is very similar to that already observed for a 2D hard disk system [22]. When linear density is low, l = 0.5, the dispersion ω_L(k) only slightly deviates from being monotonic. But as soon as linear density increases, l = 0.909, it shows well-defined minimum around wavenumber values k ~ 6 associated with position of the main peak K_max of longitudinal structure factor S(k) discussed in Figures 4, 5. The deviation from hydrodynamic dispersion law in the long-wavelength limit for both considered densities persists to be “negative.” The “negative” dispersion concerns the negative deviation of the dispersion curve ω_L(k) at the boundary of hydrodynamic regime from the linear hydrodynamic dispersion law of acoustic modes in considered quasi1-D hard disk system. These effects are similar to those observed in 2D hard disk system with density increase [22]. In contrast, the dispersion ω_T(k) of transverse excitations is essentially different from those in the case of 2D hard disk system [22]. Namely, for both linear densities l < 1 we are observing rather flat shape of the curve ω_T(k) at reduced frequency ~10 with a tendency toward higher frequency values with increase of linear density. In 2D hard disk system, the transverse excitations are of acoustic nature. Moreover, the transverse excitations are absent at low densities and there was observed a long-wavelength propagation gap when they start to appear at higher densities. It is therefore quite natural to attribute the flat transverse mode in a quasi-1D hard disk system to disks' reflections from the channel hard walls. This issue is discussed in more details below.

Parts (c) and (d) of Figure 11 show similar data for the dispersions ω_L/T(k) of longitudinal and transverse excitations but for the range of linear densities l > 1. As for the dispersions of longitudinal excitations, ω_L(k), we see the tendencies already observed in Figures 11A,B under increase of linear density, i.e., the magnitudes of ω^L(k) maxima are increasing while minima become deeper, reaching zero-frequency values at density ρ = 1.111 and being shifted toward larger wavenumber values k. The later again is consistent with the shift for the position K_max of the first peak of the longitudinal structure factor S(k) in Figures 4, 5; the “negative” dispersion in long-wavelengths region is preserved as well. Such behavior of ω_L(k) resembles one for the ordered solids and in [33] is interpreted as a consequence of the emergence a zigzag ordering in a squeezed quasi-1D system. The dispersion of transverse excitations in Figures 11C,D does not show notable changes too when the linear density was changing to l = 1.01. However, it does show dramatic changes at the highest considered linear density l = 1.111. Namely, (i) there is a sharp increase of frequency ω_T up to ~70; (ii) the dispersion curve ω_T(k) itself exhibits bubble-like shape by splitting on low- and high-frequency branches in the range of k-values that coincide with location of the maximum of dispersion ω_L(k) that implies possibility of longitudinal-transverse excitation coupling on atomic scale in a squeezed almost zigzag ordered quasi-1D hard disk system.

The MD simulations allow one to study the fluctuations of heat density in the system and their effect on collective dynamics. We would like to remind that an adiabatic propagation of sound in fluids causes small deviations of the local temperature and instantaneous temperature gradients, which give rise to relaxation processes of the local temperature via thermal diffusivity. These relaxation process is directly connected with entropy fluctuations and is responsible for the central peak of dynamic structure factors S(k, ω). In Figure 12, we show dynamic structure factors S(k, ω) as well as heat-density dynamic factors S_hh(k, ω) for three, the lowest wavenumbers in the quasi-1D hard disk system at the highest linear density l = 1.111. In both types of spectral functions, the side peaks are caused by longitudinal acoustic excitations, while the central peak for liquid state is caused by entropy fluctuations. In our case of the confined, almost zigzag structure at linear density l = 1.111 (channel length L/σ = 180 and width H/σ = 1.5), the central peaks of S(k, ω) and S_hh(k, ω) give evidence of the same temperature (sometimes called entropy) relaxation processes typical for fluid state. We calculated S_hh(k, ω) for lower densities too, and observed practically the same shape of S_hh(k, ω) but with a larger smearing of the central and side peaks in comparison with the higher densities.

In dense, nearly solid-like states of a quasi-1D hard disk system, we observed a rapid increase of the transverse frequency ω_T(k) in the long-wavelength region k ~ 0 from ω_T ~ 12 for linear density l = 1.01 in Figure 11C to the frequency ω_T ~ 70 for linear density l = 1.111 in Figure 11D. The shallow minimum in ω_T(k) profile, observed for linear density l = 1.01 at wavenumbers k ~ K_max/2, also deepens and becomes well-developed. Eventually, however, frequency ω_T(k) at k ~ K_max/2 splits into a high- and low-frequency branches at linear density l = 1.111 as shown in Figure 11D. The observed change in the dispersion ω_T(k) of transverse excitations as the system changes from rarefied to dense can be explained by formation of the zigzag structure. The transverse low frequency mode, which is due to bouncing between the two hard walls, transforms into the high frequency transverse oscillations between one wall and the nearest neighbors in the zigzag structure. This is also supported by the behavior of the transverse component of the velocity autocorrelation function ψ_yy(t) in Figure 6C. Inspired by this intriguing behavior, one of us developed an analytical theory [10] of this system, which suggests that it is related to developing window-like defects in the zigzag disks' arrangement [33]. The idea of defects of this kind (defect is a local less packing) has been introduced [12] and employed [7–9, 12–14] to describe glassy dynamics in terms of caged and uncaged states in the disk arrangements in a quasi-1D hard disk system of the channel width H/σ = 1.866. In [33], such defects were associated with the maximum contact separation of two disks along the channel, which is equal to the disk diameter σ; in [10], their distribution was found analytically as a function of the linear density. As for the densities l > 1, the horizontal contact distance between disks and the actual horizontal distance between nearest neighbor disks are very close, this distribution can be well representative of the actual disks separations Δ_x, which is supported by the computer simulation data.

At dense/close packing, the disks form perfect zigzag. As confinement weakens, the tendency to the entropy increase results in an emergence of progressively larger number of window-like defects through which pairs of next neighbor disks uncage and exchange their vertical positions (Figure 13). This theoretical prediction of [10] has been confirmed from MD simulations data for the distribution of the actual distances Δ_x between next neighbor disks (Figure 14) and is in line with the earlier result [8]. For high linear density, l = 1.111, these are distributed around Δ_x/σ ~ 0.89 close to the minimum possible contact distance along the pore, 3/2≈0.87σ. However at slightly lower linear density, l = 1.053, in addition to the widening of maximum at Δ_x/σ ~ 0.92 we see the appearance of a sharp subpeak at Δ_x/σ = 1. In contrast to the wide maximum drifting toward larger Δ_x, the subpeak always remains at Δ_x/σ = 1, although its shape is changing: it becomes more and more pronounced and finally exceeds the main peak, signaling approaching a fluid-like state. It points to the exclusive role that next neighbor disk center-to-center distance Δ_x = σ plays. Namely, creation of windows of the width of disk diameter σ (σ−windows) in the zigzag array, neither wider nor narrower, is the most effective way to gain entropy by uncaging two disks and making them extend their wondering to the whole channel width. The described mechanism of disk's caging/uncaging by their neighbors in quasi-1D hard disk system is fundamentally different to that in quasi-1D system with a long range screened electrostatic repulsion in the dusty plasma [20] where particles stay at finite distances and the zigzag can even transform into a straight line.

In a pure 1D hard disk system, the vertical motion is absent. In a densely packed quasi-1D hard disk system of the width H/σ = 1.5, it is also prevented by the full caging, but it eventually appears as the confinement weakens. As for sufficiently low density, the vertical disks' motion from one wall to another is possible, here one expects some contribution to the dispersion relation from the nearly ballistic transverse oscillation between the walls. As it comes from the maximum vertical path H − σ, this contribution ω_T1(k) to the frequency ω_T(k) at low density must be lowest possible. At high density, however, the σ−windows are rare and disks can bounce at most between one wall and the mid plane hence the lowest transverse frequency ω_T2 for higher densities is expected to be roughly twice that for low density, ω_T2 ~ 2ω_T1. The frequency ω_T2 is related to the maximal distance Δ_x at window nuclei, which require local compression and must result in high frequency longitudinal and transverse jitters. One can thus expect that the lowest transverse frequency ω_T2 and the highest longitudinal and transverse frequencies appear near same wavenumbers k. In addition, the group velocity at this wavenumber range has to be zero as windows are not transferred by the waves. The dispersion of the longitudinal and transverse excitations in quasi-1D hard disk system obtained from MD simulations are in line with this picture (Figures 11C,D). For linear density l = 1.01, the peak at Δ_x = σ (Figure 14) indicates that σ−windows in the zigzag structure are well-developed, the disk order is of a short range, and frequency ω_T1 can be identified with the practically k−independent transverse frequency ω_T(K_max/2) ~ 10 in Figure 11C. At the same time, for linear density l = 1.111, when there are almost no σ−windows (Figure 14), the transverse spectrum splits into the lowest, ω_T2 ~ 20 ≈ 2ω_T1, and the highest frequency, ω_T ~ 100, at the wavenumbers k where longitudinal frequency is maximum (see Figure 11D). At wavenumbers k ~ K_max/2, the curves ω_T1(k) and ω_T2(k) are plateaus, which indicate zero group velocities. The continuous longitudinal and transverse modes for linear density l = 1.111 are related to the short free path oscillation of mutually caged disks near the walls. Thus, the main properties of ω_T(k), which are directly related to the vertical motion, are in line with the idea of the role of σ−windows in the zigzag arrangement.

The interpretation of computer simulation data in terms of the vertical disk motion presented above is also consistent with and well-illustrated by the theoretical dependence of the total transverse force F_T(l) in Figure 2B. Indeed, as discussed above, for low linear density l, the transverse force F_T is constant implying the vertical motion to be ballistic because the free path and free time are maximum. In contrast, for high linear density l, the transverse force F_T sharply grows with l, the disks' arrangement is close to a dense zigzag, they cannot cross the middle line of the channel so that the free path and time at best are halves of their maximum values. Finally, this picture is obviously in line with the disk density distribution across the channel in Figure 3.

In conclusion, we would like to point out the novel and unexpected development, which stems from the above studies of the collective excitations. It is about a possible Kosterlitz-Thouless scenario in quasi-1D hard disk system [33]. In 2D systems, melting proceeds via the Kosterlitz-Thouless scenario [34]: a crystal develops defects, and their number is growing continuously from zero at zero temperature until the state becomes a liquid. The number of window-like defects in a zigzag arrangement behaves exactly in this way [7, 9, 10, 12]: it is zero only at dense packing and smoothly increases as density goes down. But in the Kosterlitz-Thouless scenario, it is essential that the spatial correlations decrease as a power law at high densities and exponentially at lower densities. At the same time, it is known that if the partition function of a system possesses a transfer matrix property, which has been widely accepted to be the case for a quasi-1D hard disk system, then the correlations can only decay exponentially. Our search for a different correlation behavior is motivated as follows. First, our already published [33, 35] and tentative molecular dynamic results on the system described here indicate that a power law decay is possible. Second, it is shown in Appendix in [10] that the transfer matrix property may be not so universal for quasi-1D systems, which leaves a room for alternative theory. The work is in progress.

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

AH and AT designed the model and the computational framework. AH performed the event driven MD simulations to generate particle trajectories. AT performed the event driven MD simulations of particle distributions. TB used MD-generated trajectories and performed the numerical calculations of collective dynamics properties. VP has contributed the interpretation of the transverse excitation modes in terms of the windowlike defects and the connection between the transverse disk motion and the total forces and pressures on the channel's boundaries. AT supervised the project and wrote the manuscript with support from VP, TB, and AH. All authors discussed the results and contributed to the final manuscript.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.