Static and Dynamic Properties of Two Dimensional Coulomb Clusters

TL;DR

This study investigates how two-dimensional Coulomb clusters transition from solid-like to liquid-like states with temperature, revealing glassy dynamics and long-lived inhomogeneities influenced by confinement irregularities.

Contribution

It provides new insights into the static and dynamic behaviors of Coulomb clusters, highlighting glassy signatures and long-range correlations in irregular traps.

Findings

Static correlations show hexatic-like phase at low temperature

Particle dynamics slow down, resembling supercooled liquids

Long-lived inhomogeneities due to confinement and Coulomb interactions

Abstract

We study the temperature dependence of static and dynamic responses of Coulomb interacting particles in two-dimensional traps across the thermal crossover from an amorphous solid- to liquid-like behaviors. While static correlations, that investigate the translational and bond orientational order in the confinements, show the footprints of hexatic-like phase at low temperature, dynamics of the particles slow down considerably in this state -- reminiscent of a supercooled liquid. Using density correlations, we probe intriguing signatures of long-lived inhomogeneities due to the interplay of the irregularity in the confinement and long-range Coulomb interactions. The relaxation at multiple time scales show stretched-exponential decay of spatial correlations in irregular traps. Temperature dependence of characteristic time scales, depicting the structural relaxation of the system, show striking similarities with those observed for the glassy systems indicating that, some of the key signatures of supercooled liquids emerge in confinements with lower spatial symmetries.

Tables

Table 1. Table 1: Estimation of the parameter Γ = π ​ n / T X Γ 𝜋 𝑛 subscript 𝑇 𝑋 \Gamma=\sqrt{\pi n}/T_{X} , characterizes the melting in bulk 2 ​ D 2 D \mathrm{2D} systems, for irregular and circular confinements. For bulk systems Γ ≈ 137 Γ 137 \Gamma\approx 137 at the transition. Note that the average particle density n 𝑛 n is estimated by computing the area covered by N 𝑁 N particles in the respective confinements.

$N$	Irregular			Circular
	$n$	$T_X$	$\Gamma$	$n$	$T_X$	$\Gamma$
75	2.42	0.017 $\pm$ 0.003	162	1.90	0.011 $\pm$ 0.001	222
150	2.53	0.019 $\pm$ 0.001	148	2.12	0.019 $\pm$ 0.002	136
500	3.01	0.021 $\pm$ 0.001	146	2.78	0.020 $\pm$ 0.001	148

Table 2. Table 1: Details of the fitting parameters, τ g subscript 𝜏 𝑔 \tau_{g} and β 𝛽 \beta , of the cage correlation function, c g ​ ( t ) subscript 𝑐 𝑔 𝑡 c_{g}(t) , for circular and irregular confinements. For both the confinements the decay of c g ​ ( t ) subscript 𝑐 𝑔 𝑡 c_{g}(t) is stretched exponential.

Temperature	Irregular		Circular
	${\tau }_g$	$\beta$	${\tau }_g$	$\beta$
0.010	2271.83	0.446	783.26	0.508
0.020	188.002	0.442	247.56	0.574
0.030	67.6891	0.553	34.035	0.639
0.050	19.8172	0.596	15.638	0.603

Equations

$$\mathcal{H}=\frac{q^2}{ϵ}\sum _{i<j=1}^N\frac{1}{|{\vec{r}}_i-{\vec{r}}_j|}+\sum _{i=1}^N{V}_{\mathrm{conf}}(r_i),$$

$$V_{\mathrm{conf}}^{\mathrm{Cr}}(r)=\alpha r^2,\text{ where}\alpha =m{\omega }_0^2/2,$$

$$V_{\mathrm{conf}}^{\mathrm{Ir}}(r)=a\{x^4/b+b{y}^4-2\lambda x^2{y}^2+\gamma (x-y)xyr\}.$$

$${\mathcal{H}}^{\mathrm{Cr}}=\sum _{i<j=1}^N\frac{1}{|{\vec{r}}_i-{\vec{r}}_j|}+\sum _i{r}_i^2$$

$${\mathcal{H}}^{\mathrm{Ir}}$$

$${\psi }_6(k)=\frac{1}{N_b}\sum _{l=1}^{N_b}{\mathrm{e}}^{i6{\theta }_{kl}}$$

$${\varphi }_6(k)=\left|{\psi }_6(k{)}^{\ast }\frac{1}{N_b}\sum _{l=1}^{N_b}{\psi }_6(l)\right|$$

$${\chi }_{\psi }$$

$${\chi }_{\varphi }$$

$$g_6(t)=⟨{\psi }_6^{\ast }(t){\psi }_6(0)⟩$$

$$Q(t)=⟨\frac{1}{N}\sum _{i=1}^{N}w(|{\vec{r}}_i(t_0+t)-{\vec{r}}_i(t_0)|)⟩$$

$${\chi }_4(t)=\frac{1}{N}[⟨Q^2(t)⟩-⟨Q(t){⟩}^2]$$

$${\mathbf{L}}_{\mathbf{i}}(\mathbf{t})=f(r_{ij})$$

$$C_g(t)=\frac{⟨{\mathbf{L}}_{\mathbf{i}}(\mathbf{t})\cdot {\mathbf{L}}_{\mathbf{i}}(0)⟩}{⟨{\mathbf{L}}_{\mathbf{i}}^2(0)⟩}$$

$$\frac{{\sigma }_i(t)}{\Delta r_{ij}(t)}\le d;\text{ for at least}j=3$$

$$c_g(t)\propto \exp \left(-{\left(\frac{t}{{\tau }_g}\right)}^{\beta }\right)$$

Full text

I INTRODUCTION

Long-range interacting classical particles in disordered media can give rise to exotic phases, e.g, hexatic glass Chudnovsky (1989), which is characterized by short range positional order and quasi-long range orientational order. This is quite similar to the usual hexatic phase in two dimension ($\mathrm{2D}$), as discussed in Kosterlitz, Thouless Kosterlitz and Thouless (1972, 1973), Halperin, Nelson, and Young (KTHNY) Halperin and Nelson (1978); Nelson and Halperin (1979); Young (1979) melting scenario. Though a true hexatic glass in the thermodynamic limit attracted critical debates Toner (1991a); Chudnovsky (1991), it may be realized in finite systems Toner (1991b). However, there have been attempts made to understand experiments with colloidal particles Kusner et al. (1994), binary-alloy Zippelius et al. (1980), magnetic bubble lattice Seshadri and Westervelt (1991) and disordered type II superconductors Chudnovsky (1989) in light of a hexatic glass.

Confined systems, with small number of particles, are of great current interests due to their controlled experimental tunability over a wide range of parameters. Such systems are not only significant from technological perspective, they are also very important for fundamental Physics: They are ideal play ground for exploring complex interplay of interactions and disorder. In this paper, we develop an understanding of such interplay by studying static and dynamic correlations in nano-clusters with (long-range) Coulomb-interacting particles.

Standard order-disorder transitions do not occur in $\mathrm{2D}$ Mermin and Wagner (1966) because, a true ordered phase is weakened by fluctuations producing only quasi-ordered state. Still a transition from solid to liquid occurs with an intervening hexatic phase which is characterized by short range positional order and quasi-long range bond-orientational order. Such two step melting scenario, predicted by the KTHNY theory for clean $\mathrm{2D}$ systems, is realized in certain experiments Han et al. (2008); von Grünberg et al. (2004); Peng et al. (2010) and simulations Prestipino et al. (2011); Kapfer and Krauth (2015). The possibility of extending such ideas in the presence of disorder is still a subject of current research Nelson (1983); Serota (1986); Deutschländer et al. (2013); Nattermann et al. (1991). Does the melting of Coulomb interacting particles in $\mathrm{2D}$ system corroborate with the KTHNY melting, even in the absence of any disorder? In spite of thirty years of research, there is no definitive answer, to the best of our knowledge Mazars (2015).

Study of dynamics across the melting often provides crucial additional insights Zahn and Maret (2000). The dynamical behavior of a system indeed helps in probing the intricacies of different thermal phases. The dynamical response of $\mathrm{2D}$ systems, close to the liquid-hexatic transition, show some striking similarities with that of liquids close to the glass transition Ediger (2000). Molecular dynamics simulation of $\mathrm{2D}$ system of colloid particles Zangi and Rice (2004) and experiment with granular materials Dauchot et al. (2005); Reis et al. (2007) show that with the increase in orientational order, the dynamics of the particles slow down and become heterogeneous. Particles undergo ‘caging effect’ (i.e., getting temporarily trapped in the cage of its neighbors) and cooperative dynamics. These results into sub-diffusive behavior of the mean square displacements as observed in experiments Weeks et al. (2000) and simulations Kob and Andersen (1995a, b); Donati et al. (1998).

While true phase transitions pertain to bulk systems only, for finite systems the notion of solid- and liquid-like ‘phases’ have been used successfully to characterize qualitative behaviors of macroscopic systems Böning et al. (2008). In recent years, static and dynamic properties of $\mathrm{2D}$ finite clusters were studied for different types of confinements Bedanov and Peeters (1994); Bhattacharya and Ghosal (2013); Bhattacharya et al. (2016); Ash et al. (2016); Kong et al. (2003) and interaction potentials Apolinario et al. (2005); Bhattacharya and Ghosal (2013); Girotto et al. (2014) across the thermal melting. These theoretical studies are motivated by experimental realizations, such as colloidal suspensions Bubeck et al. (1999); Costa Campos et al. (2013); Golosovsky et al. (2002), confined plasma Melzer et al. (2012); Thomsen et al. (2012); Chu and I (1994), electrons in quantum dots in high magnetic fields Ashoori (1996); Kouwenhoven et al. (2001), radio-frequency ion traps Buluta and Hasegawa (2009) and electrons on the surface of liquid helium Grimes and Adams (1979). Recent studies of static Bhattacharya and Ghosal (2013) and dynamic Ash et al. (2016) properties of Coulomb interacting particles in irregular geometry have confirmed the thermal crossover from ‘solid’ to ‘liquid’-like phases. From the study of spatio-temporal correlations Ash et al. (2016) in such system, it is found that with decrease in $T$ some of the key signatures of glassy dynamics emerges. For example, as system approaches solidity from liquid-like phase, the dynamics of the particles become slow and heterogeneous. Is this a signature of a hexatic-glass phase?

We address this broad question by presenting in this paper the results of extensive computer simulations on the static and dynamic properties of Coulomb interacting particles in different confinement geometries. From our analysis of static properties, we find that while the positional order is depleted even at the lowest $T$ for an irregular trap, a solid-like phase can still be identified at low $T$ due to the bond orientational order. A liquid-like phase emerges Bhattacharya and Ghosal (2013) with the breaking down of bond orientational order beyond a crossover temperature $T_{X}$, however, the inhomogeneities in the liquid persists up to a much larger temperature – about an order of magnitude larger than $T_{X}$. Beyond this large temperature scale ($\sim 10T_{X}$), the system crosses over to a standard isotropic liquid. Our results are indicative of a hexatic-like low $T$ phase in the confinements. Finally, the analysis of the trajectories of individual particles and the $T$-dependence of relaxation time-scales show characteristics quite similar to those observed for glass formers near the glass transition temperature Ediger (2000).

The rest of the paper is organized as follows: In Sec. II we will give the details of the models and methods used in our simulations emphasizing the way disorder is introduced in finite systems. We will also discuss about the methods used to study the system. Then we present, in Sec. III, our results for static properties where we analyze the thermal crossover through the temperature dependence of bond orientational order. In Sec. IV dynamic properties are analyzed to identify, in particular, any signature of a glass-like behavior, and also to study the temperature dependence of structural relaxation time(s). Finally, we conclude in Sec. V.

II MODEL and Method

We consider $N$ classical particles each with charge $q$ in a confining potential $V_{\rm conf}(r)$. These particles interact via long range Coulomb potential and are restricted to move in $\mathrm{2D}$. The Hamiltonian describing such a system reads as:

[TABLE]

where, $r_{i}=|\vec{r}_{i}|=\sqrt{x_{i}^{2}+y_{i}^{2}}$ is the distance of the $i$-th particle from the origin. Here, first term in the Hamiltonian represents Coulomb repulsion between particles in a medium with dielectric constant $\epsilon$. We consider long range Coulomb repulsion alone, because systems with small number of particles, particularly, in the presence of disorder, are expected to offer a very weak screening, if at all. In our study we considered two types of confinement potential; (a) parabolic (having circular symmetry):

[TABLE]

and (b) Irregular (that lacks all spatial symmetries) Bohigas et al. (1993); Bhattacharya and Ghosal (2013):

[TABLE]

The solid phase of Coulomb clusters at low $T$ is called a Wigner Molecule (WM) Egger et al. (1999), because it mimics the physics of Wigner crystal Wigner (1934). In finite systems, we refer a WM as circular Wigner molecule (CWM) in a circular trap, and as irregular Wigner molecule (IWM) in an irregular trap.

We rescale the length $r^{\prime}\rightarrow\phi^{1/3}\alpha^{-1/3}r$ and energy $E^{\prime}\rightarrow\phi^{2/3}\alpha^{1/3}E$, where $\phi=q^{2}/\epsilon$, in such a way that the CWM Hamiltonian transforms to Bedanov and Peeters (1994):

[TABLE]

correspondingly the time scale will also be renormalized to $t^{\prime}=\hbar~{}\phi^{-2/3}\alpha^{-1/3}t$. Such renormalization of length, energy and time will eventually make temperature to be expressed as $T^{\prime}=E^{\prime}/k_{B}$, where $k_{B}$ is the Boltzmann constant. In order to have an estimate of these new length and time scales in conventional units, we consider electrons in GaAs heterostructure with a typical confinement energy of $\hbar\omega_{0}=1$ meV. The above scaled length $(r^{\prime})$, energy $(E^{\prime})$ and time $(t^{\prime})$ unit takes the value of $630$ Å, $1.7$ meV and $376$ fs, respectively. In extracting these values we have used, $q=-e$, the charge of an electron, $m=0.067m_{e}$, the mass of an electron and $\epsilon=13$ Reed and Kirk (1989).

$V_{\rm conf}^{\rm Cr}$ is quadratic in length-scale while $V_{\rm conf}^{\rm Ir}$ is quartic. In order to facilitate a justified comparison between the two, we need to express the quartic irregular confinement in the units of quadratic circular confinement. This is achieved by setting the parameter $a=\left(m\omega_{0}^{2}/2r^{\prime 2}\right)a^{\prime}$, which includes both $\omega_{0}$ and a scaling factor $a^{\prime}$, for irregular confinement. Thus a judicious choice of $a^{\prime}$, which now controls the strength of the irregular confinement, brings the two confinements on equal footing on dimensional ground. Following the same scale transformations as for the circular case, we obtain the following Hamiltonian for $V_{\rm conf}^{\rm Ir}$:

[TABLE]

There is still a parameter, $a^{\prime}$, left to be fixed for $V_{\rm conf}^{\rm Ir}(r)$. For a given $N$, we choose the value of $a^{\prime}$ according to the following prescription: We know that the model of a $\mathrm{2D}$ electron gas, neutralized by a uniform positive background, and interacting via Coulomb interaction can be characterized by a single dimensionless quantity $\Gamma=\sqrt{\pi n}/T$, where $n$ is the average number density Gann et al. (1979). Thus, for a given $T$, the thermodynamic properties of the system is completely determined by its density, $n$. In the present context of finite systems, we adopt the same view and assume that our systems are also characterized by the same parameter $\Gamma$ (which we will justify below). In order to make a meaningful comparison between systems in different confinements but at same $T$ and with same $N$, we fix the density by tuning $a^{\prime}$ for irregular confinement and $\alpha$ for circular confinement. The assumption of single parameter description in terms of $\Gamma$ is later indicated (in sec. IIIC) by its bulk value at melting, irrespective of $N$ and trap geometries. Note that $a^{\prime}$ tunes the average density by making the quartic oscillator narrow or shallow. For a given $N$, we fixed the value of $a^{\prime}$ for the irregular confinement in such a way that the average inter-particle distance is equal to that for circular confinement at the lowest $T$. In this paper, we report results for $N=75,150$ and $500$ number of particles which require $a^{\prime}=0.10,0.055$ and $0.020$, respectively, following above procedure.

Note that while the thermal crossover in circular trap is already studied in great details Bedanov and Peeters (1994), the main objective of this paper lies with the understanding of melting in irregular confinement. We study circular confinement mainly for comparison and thus move those results predominantly in supplementary materials.

The parameter $\lambda$ in irregular confinement controls chaoticity in the single-particle dynamics Bohigas et al. (1993). Tuning $\lambda$ from zero to unity generates periodic to chaotic motion of a sole particle in the trap. Chaotic motion, along with broken spatial symmetries are taken as the signatures of disorder in our study. The parameter $b(=\pi/4)$ breaks the symmetry of a square and $\gamma$ breaks the reflection symmetry. We consider $\lambda\in[0.565,0.635]$ and $\gamma\in[0.10,0.20]$ Jiang et al. (2003). The values of different parameters are adjusted to generate self-similar copies of the system over which statistics are collected for the purpose of “disorder averaging” of physical observables Ullmo et al. (2003). Results from CWM are averaged for many independent simulations.

To study the static properties, we carried out (classical) Metropolis Monte Carlo (MC) simulation Frenkel and Smit (2001) aided by simulated annealing algorithm Kirkpatrick et al. (1983). Dynamical responses are studied using molecular dynamics simulation (MD) Frenkel and Smit (2001). To achieve a desired $T$, we have used Berendsen like thermostat  Frenkel and Smit (2001) during the equilibration. Once equilibration is done, we have used conventional velocity-Verlet algorithm Frenkel and Smit (2001) to integrate the equations of motion. We have performed MD runs up to $2\times 10^{6}$ steps with a time step size of $dt=0.005t^{\prime}$. We find excellent match of different observables obtained from Monte Carlo simulation with the time-averaged results of those quantities obtained from molecular dynamics simulation, except for the lowest temperatures featuring slow dynamics, as we discuss below. Such match of the physical results from the two independent methods of numerical calculations validate the correctness of our findings.

III STATIC PROPERTIES

The broken symmetry state of a system at low $T$ can be identified by its order. For example, crystalline solid is characterized by the long range positional and orientational order. We thus proceed to explore first the positional order in our confined systems at the lowest $T$. Note that even though translational symmetry is broken by confinements, a circularly confined system would still possess an azimuthal periodicity Bedanov and Peeters (1994).

III.1 POSITIONAL ORDER

In order to quantify the positional order in IWM, we plot the Fourier component of local density, $\rho(\vec{k})=\sum_{j=1}^{N}\mathrm{exp}[i\vec{k}\cdot\vec{r}_{j}]$, in Fig. 1(b). Here, $\vec{k}$ is the momentum vector and $\vec{r}_{j}$ represents the position vector of the $j$-th particle. In all analysis, distances are measured in the unit of $r_{0}$, the mean inter-particle spacing, at the lowest $T$, between neighboring particles (which is same by construction for irregular and circular confinements, for a given $N$). For a perfect crystalline structure, $\rho(\vec{k})$ consists of Bragg peaks at the reciprocal lattice vectors and such peaks broaden with $T$, finally disappearing upon melting. Fig. 1(a) shows the ground state configuration, obtained following the prescription of Ref.Bhattacharya and Ghosal (2013), of $N=500$ particles in irregular confinement (with $\lambda=0.635,\gamma=0.20$) and Fig. 1(b) shows $\rho(\vec{k})$ corresponding to this configuration. We find no strong peak in reciprocal space, apart from the one at $\vec{k}=0$. Within the diffuse pattern at $\vec{k}\neq 0$, we broadly identify six patches. Broadness of these patches ensures that there is, at best, a very short range (much shorter than the system size) positional order, if any, in IWM. We can probe this order further by considering the fact that the ground state configuration for a trapped system is a result of the interplay between the formation of a triangular lattice (minimum energy configuration for a bulk $\mathrm{2D}$ Coulomb system Bonsall and Maradudin (1977)) towards the central part of the trap (i.e., far away from the boundary) and the geometry of the confining potential near the boundary. In search of any positional order deep inside the system, we consider only those particles which belong to a small region (of length $l_{x}=7r_{0}$ along $x$-axis and $\sqrt{3}l_{x}/2$ along $y$-axis, to accommodate a commensurate triangular lattice structure) near the center of the confinement. These particles are colored orange in Fig. 1(a). We then calculate $\rho(\vec{k})$ considering only the orange particles. Fig. 1(c) shows that the six patches become somewhat prominent enhancing positional order in this subsystem. This is also obvious from Fig. 1(a), where the lattice lines become wiggly even in the subsystem due to irregular geometry. At $T=0.020$, which characterizes the crossover temperature (as discussed later), $\rho(\vec{k})$ shows a diffuse pattern ( Fig. 1(d) ), similar to that at $T=0$, establishing little evolution of positional order, if any, with temperature. Similar analysis for the CWM shows stronger positional order, primarily reflecting the azimuthal symmetry, though it gets depleted upon averaging over whole system (see Fig. 1 in supplementary materials).

In the absence of any appreciable positional order, any solidity in IWM is contributed by the orientational order in 2D. Motivated by this assertion, we next analyze the bond orientational order in the system.

III.2 BOND ORIENTATIONAL ORDER (BOO)

In Fig. 1(a) we see that particles which are not on the boundary are mostly surrounded by six nearest neighbors oriented in a hexagonal fashion and thus indicate six-fold bond orientational order (BOO) of a triangular lattice. In $\mathrm{2D}$ systems, the BOO parameter, $\psi_{6}(k)$, for $k$-th particle is defined as Nelson (2002):

[TABLE]

Here, $\theta_{kl}$ is the angle made by the bond between the particles $k$ and $l$ with an arbitrary axis and $N_{b}$ is the number of nearest neighbors of particle $k$. We identify the nearest neighbours of a given particle by Voronoi construction Tipper (1991). $\left|\psi_{6}\right|$ achieves a maximum value of unity in a perfect triangular lattice and would get smaller with distortions.

Fig. 2, shows the $T$-dependence of the probability distribution, $P(\left|\psi_{6}\right|)$, of $\left|\psi_{6}\right|(=\left|\psi_{6}(k)\right|$ for $k=1,2,\cdots N)$ for IWM with $N=500$ (panel a), $150$ (panel c) and $75$ (panel e) particles. For $N=500$ (Fig. 2(a)), this distribution is sharply peaked around unity at low $T(=0.002)$, describing an orientationally ordered phase at low $T$. With increase in $T$, the peak broadens due to thermal fluctuations that distorts the bond angles $(\theta_{kl})$ with the neighbors. For $T>0.02$, we get a very broad distribution signifying the breakdown of orientational order. Thus, the $T$-dependence of $P(\left|\psi_{6}\right|)$ signals a thermal crossover from an orientationally ordered solid-like to a disordered liquid-like phase in confined systems. It is interesting to note that $P(\left|\psi_{6}\right|)$ shows hardly any thermal evolution for $T>0.03$. This is found true irrespective of the confinement geometry or total number of particles,$N$. This insensitivity confirms that the scrambling of orienational order and hence melting is complete by $T\approx 0.03$.

With decreasing $N$ (Fig. 2(c) for $N=150$ and Fig. 2(e) for $N=75$), we notice a strengthening of a small and feeble peak in $P(\left|\psi_{6}\right|)$ near $\left|\psi_{6}\right|=1$. We assert that such peak is spurious. We verified that this is contributed predominantly by the particles on the boundary, having very few nearest neighbors. Note that a particle with a sole neighbor must yield $\left|\psi_{6}\right|=1$, irrespective of the value of bond angle. Decrease in $N$ increases the ratio of the number of boundary- to bulk-particles, enhancing the strength of such spurious peak. From our overall study, we find that systems with $N\leq 100$ tend to show non-universal behavior, not only for $P(\left|\psi_{6}\right|)$ but also for other observables we describe below. On the other hand, systems with $N=100$-$500$ describe generic behavior featuring qualitatively similar results.

The nature of melting by disrupting bond orientational order is usually tracked by studying the large $r$ behavior of bond-orientational correlation function, $g_{6}(r)=\langle\psi_{6}(r)\psi_{6}^{*}(0)\rangle$. However, for finite systems with small number of particles, $g_{6}(r)$ is not very useful, because any trend in its $r$-dependence is hard to discern Bhattacharya and Ghosal (2013) due to limited values of $r$.

III.3 PROJECTION OF BOND ORIENTATIONAL ORDER

In order to get additional insight into the degree of local orientational order, we, next, look into the short distance projection of bond-orientational correlation function. The magnitude of the projection of bond orientational order, $\phi_{6}(k)$, of particle $k$ onto the average local bond orientational order of its nearest neighbours is defined as Larsen and Grier (1996); Dillmann et al. (2012):

[TABLE]

While $\left|\psi_{6}(k)\right|$ measures how the neighbors of particle $k$ adjust locally on a hexagonal environment, $\phi_{6}(k)$ determines how the orientation of a particle fits into the mean orientational field of its nearest neighbors. Thus, $\phi_{6}(=\phi_{6}(k)$ for $k=1,2,\cdots N)$ represents the bond orientational correlation up to next nearest neighbor distances.

At low $T$, as most of the particles are surrounded by six regular neighbors, and we expect $\phi_{6}\sim 1$. With increasing $T$, due to depleted orientational order, $\phi_{6}$ decreases and eventually becomes close to zero at high $T$. Thus, $\phi_{6}$ distinguishes a high $T$ phase clearly from a low $T$ phase.

In Fig. 2(b,d and f) we show the distribution of $\phi_{6}$ at different $T$ for IWM with $N=500,150$ and $75$, respectively. At low $T(=0.002)$, the distribution is sharply peaked at $\phi_{6}\approx 1$ (see Fig. 2(b)), depicting a strong correlation of the orientational ordering of a particle with its nearest neighbours. As $T$ increases, such correlation reduces and at high $T(=0.050)$, the distribution becomes sharply peaked around $\phi_{6}\sim 1$. Thus, for $T>0.02$, the local six-fold order is lost even at the second shell, representing a ‘liquid’-like state. We see that with decrease in $N$ (Fig. 2(d,f)), the sharpness of the peak around $\phi_{6}\sim 1$ gets considerably reduced at the lowest $T(=0.002)$ as the effect of the boundary particles becomes more prominent for smaller system size. For $N=75$ (Fig. 2(f)), an orientational ordering is hard to discern up to the second neighbour even at low $T$.

An earlier investigation Bhattacharya and Ghosal (2013) of IWM with $N=150$ found $T_{X}\approx 0.02$ (for our parameters) by studying the $T$-dependence of specific heat and Lindemann ratio, among other quantities. Here we identify $T_{X}$ by studying the fluctuation of the orientational order parameters $|\psi_{6}|$ and $\phi_{6}$, defined as

[TABLE]

Here, $\chi$’s are the generalized susceptibility correspond to order parameter $\left|\psi_{6}\right|$ and $\phi_{6}$, respectively. In Fig.3, we show the temperature dependence of $\chi_{\psi}$ (panel a,b) and $\chi_{\phi}$ (panel c,d) for different system sizes. While, for a bulk system, susceptibility diverges at the critical point, such divergence turns into a crossover in a finite system. The fluctuations in $\chi_{\psi}$ and $\chi_{\phi}$ attain their maxima at the crossover temperature $T_{X}$ distinguishing the low temperature solid- and high temperature liquid-like phases, as seen from Fig. 3(a-d). The peak at $T_{X}\sim 0.02$ becomes sharper with increasing $N$ (though for $N=75$, the circular confinement shows deviation), and we generally find the value of $T_{X}$ quite insensitive to the nature of confinement. Thus, $T_{X}$ estimated from bond orientational order parameter is in good agreement with its previous estimates Bhattacharya and Ghosal (2013).

In order to compare $T_{X}$ with the corresponding melting temperature for bulk $\mathrm{2D}$ systems, we estimate the parameter $\Gamma(=\sqrt{\pi n}/T)$ at $T=T_{X}$. It is found that for bulk systems, melting occurs for $\Gamma\approx 137$ Gann et al. (1979). In Table 1, we show $\Gamma$ for different $N$ and found that $\Gamma$ remains close to $137$ within acceptable tolerance for the large $N$’s of our study.

Since $\phi_{6}$ is a projection of $\psi_{6}$, $\phi_{6}(k)\leq\left|\psi_{6}(k)\right|$. Because maximum value of either of $\psi_{6}$ or $\phi_{6}$ is unity, $\phi_{6}(k)+\left|\psi_{6}(k)\right|\leq 2$. Hence, we can estimate melting by describing each particle to be either ‘solid’- or ‘liquid’-like using the majority rule Larsen and Grier (1996): particle $k$ is solid-like if $\phi_{6}(k)+\left|\psi_{6}(k)\right|>1$ and liquid-like otherwise. In Fig. 3(e-f) we have shown the fraction of solid- and liquid-like particles at different $T$ for irregular and circular confinements, respectively. We see that at very low $T(=0.002)$ almost $90\%$ of the particles are solid like while at high $T(=0.050)$ particles are mostly liquid like. Interestingly, the two curves cross each other around the same $T_{X}$ as obtained from other quantities above.

Few points deserve special mention here: A phase with depleted positional order but featuring an orientational order in the scale of system size is indicative of a hexatic phase within the KTHNY description of bulk 2D melting. However, in systems with disorder, like the irregular confinement considered here, similar characteristics should identify a so called ‘hexatic glass’ phase. This phase differs from regular hexatic phase in the sense that it supports both the thermal and quenched defects, whereas the hexatic phase in KTHNY theory allows only thermal defects. While the existence of a thermodynamic hexatic glass phase has been debated Toner (1991a), its signature in small finite systems, like IWM, can possibly be discerned as our results in preceding section suggest.

IV Dynamic properties

A slow and heterogeneous motion of constituent particles is a hallmark of glassy systems. Below, we resort to dynamical characterization of different phases of our system in order to explore its glassy nature, if at all. Some aspects of intriguing slow and heterogeneous dynamics in IWM, primarily extracted from the analysis of Van Hove correlation function, are recently reported Ash et al. (2016) by us. Here we expand on those results, firstly, by characterizing qualitatively the motional signatures in the traps, and subsequently, by extracting time scales associated with the structural relaxation causing the depletion of bond orientational order.

IV.1 Analysis of the trajectory

Focusing on a specific realization of the irregular confinement, $\{\lambda,\gamma\}=\{0.635,0.20\}$, we show in Fig. 4(a) the initial configuration of $N=150$ particles (thick dots) after MD equilibration. On this we also superimpose the trajectory of three particles (indexed as $1,2$ and $3$) during the whole MD-simulation. We see that, at low $T$, the nature of the dynamics of individual particles can be broadly classified in following two categories:

(i) There are particles, like the one indexed $1$ in Fig. 4(a), which only rattles around their equilibrium position over the duration of whole MD simulation.

(ii) Additionally, there are particles, indexed $2$ and $3$ in Fig. 4(a), which move longer distances (by several $r_{0}$) in the course of time. Their trajectories are enlarged in Fig. 4(b) and (c) (along with a solid line representing the displacement). These particles jiggle around their equilibrium position for certain time and then jump close to location of its neighbor and spend some time before jumping to previous (reversible motion) or a new location (irreversible motion). Such occasional jumps can be captured by looking into the instantaneous position, $r_{i}(t)$ in Fig. 4(d) where $i$ represent the particle index. We see that within a sufficiently small time window these particles change their positions moving by a distance of the order of $r_{0}$. This motion can be regarded as the cage breaking process and avalanches of such cage-breaking events give rise to string-like pattern for the displacement $\Delta\vec{r}_{i}(t)=\{\vec{r}_{i}(t)-\vec{r}_{i}(0)\}$ of the particles as shown in Fig. 4(e-f) for $N=500$. We, however, note that the length of such string in the scale of system size, decreases by going from $N=150$ to $N=500$. This cooperative motion in irregularly trapped Coulomb clusters is similar to those observed in other systems such as colloids Nagamanasa et al. (2011), Lennard-Jones mixture Donati et al. (1998) and $\mathrm{2D}$ dusty plasma Chan et al. (2009). Equivalent analysis for circular confinement (Fig. 3 in supplementary materials) also shows heterogeneous dynamics for the particles but the spatial heterogeneity in that case is dictated by its azimuthal symmetry. We conclude this section by asserting that similar motional footprints have been found for each realization (identified by a specific pair of $\{\lambda,\gamma\}$) of irregular confinement we studied.

IV.2 Inhomogeneity of the liquid beyond the crossover

In a recent study Ash et al. (2016), we found that the heterogeneous nature of the dynamics becomes profound near the thermal crossover, where the self-part of Van Hove correlation $(G_{s}(r,t))$ Hansen and McDonald (2006) shows stretched exponential spatial decay ($G_{s}(r,t)\sim\mathrm{exp}(-lr^{k})$). Here, we extend such studies to track the $N$-dependence of $k(t)$ for different $T$ as shown in Fig. 5 (panel a: $N=500$, panel b: $N=150$ and panel c: $N=75$) for $V_{\rm conf}^{\rm Ir}$. Here, we denote $k$ for $V_{\rm conf}^{\rm Ir}$ and $V_{\rm conf}^{\rm Cr}$ as $k_{\rm Ir}$ and $k_{\rm Cr}$, respectively. While the stretched exponential behavior is generic for all $N$, the finer details differ. In particular, the recovery of the final isotropic liquid like behavior signaled by $k_{\rm Ir}\approx 2$ for large $t$ occurs differently for different $N$. In Fig. 5(c), we see that an isotropic liquid is not at all recovered for $N=75$ up to the largest time of our simulation, even at high $T(=0.10)$. System with $N=500$ particles (Fig. 5(a)) attains such limit more smoothly than $N=150$ (Fig. 5(b)), indicating a surprisingly large temperature window of inhomogeneous liquid (as signaled by $k_{\rm Ir}<2$) will gradually shrink as our nano-cluster approaches the bulk limit. We find that $k_{\rm Ir}$ attains the minimum value, which increases with $N$, at $T\sim T_{X}$. Our results for low $T$ indicates that the conjecture of universal exponential tails of $G_{s}(r,t)$ in glass formers Chaudhuri et al. (2007) hold even in case of Coulomb clusters.

Similar analysis for circular confinement shows that while $k_{\rm Cr}\approx 1$ at low $T$ , $1<k_{\rm Cr}<2$ at intermediate $T$, for all $N$, implying a stretched Gaussian decay of $G_{s}(r,t)$ in contrast to stretched exponential decay observed in irregular confinement (for details see Fig. 4 in supplementary materials). So, we see while the ‘solid’ phase is quite similar in circular and irregular geometry, the Coulomb liquid is quite inhomogeneous in irregular geometry for a wide range of $T(0.02<T<0.25)$.

IV.3 Temporal bond orientation correlation function

Recent experimental advances Deutschländer et al. (2013); Zahn and Maret (2000) show that study of long time dynamical behavior of a system at different $T$ can yield crucial informations in identifying the solid, liquid and even hexatic phases for $2D$ systems. Here, we analyze bond orientational correlation function in time Deutschländer et al. (2013), defined as

[TABLE]

In the solid phase $g_{6}(t)$ remains constant in $t$ while it decays exponentially in the isotropic liquid phase for bulk systems. In the hexatic phase, $g_{6}(t)$ shows an algebraic decay, i.e., $g_{6}(t)\sim t^{-\alpha(T)}$ and $\alpha(T)\sim 1/8$ Nelson (2002); Zahn and Maret (2000) at the hexatic to liquid transition, though the transition turns into a crossover with disorder Nattermann et al. (1991).

Fig. 6 shows the $t$ dependence of $g_{6}(t)$ at different $T$ with $N=500$ for IWM (panel a) and CWM (panel b) which clearly identifying three distinct behavior of $g_{6}(t)$ as mentioned: At very low $T$ ($\leq 0.006$), $g_{6}(t)$ remains nearly flat in $t$ implying a ‘solid-’ like phase. It decays exponentially in the high $T(\geq 0.030)$ liquid phase whereas in the intermediate $T$ the decay is algebraic $(g_{6}(t)\sim t^{-\alpha(T)})$. $\alpha(T)$ shows a very weak $T$ dependence in irregular confinement and $\alpha(T_{X}\simeq 0.02)=0.17$, a value which is larger than the KTHNY prediction. On the contrary $\alpha(T)$ shows a stronger $T$ dependence in circular confinement with $\alpha(T_{X})=0.15$. Such behavior of $g_{6}(t)$, qualitatively discerning its power law decay, is a compelling evidence for the existence of a hexatic glass like phase in our irregularly confined Coulomb clusters. Similar inference holds for $N=150$ and is illustrated in Fig. 5 in supplementary materials.

IV.4 Time scales of structural relaxation

Signatures of multi-scale temporal relaxations in confined systems is recently elucidated Ash et al. (2016) in our earlier study. Here, we focus on extracting relevant time scales for structural relaxation from the overlap function, cage correlation function and persistence and exchange times. We will discuss how these signify the underlying motional processes quantifying the nature of the dynamics. We will finally look into the cross-correlation between these times scales to understand the inter-relations between these time-scales.

IV.4.1 Structural Relaxation Time from Overlap Function

We estimate the time evolution of the self part of a two-point density correlation function, called the overlap function $Q(t)$ Karmakar et al. (2014, 2009), defined as:

[TABLE]

where $w(r)=1.0$ if $r<r_{c}$ and zero otherwise. The angular parentheses denotes averaging of results over the time origin, $t_{0}$, and also over different realizations of the disorder. Since particles mostly rattle in a small distance around their equilibrium positions (within $r_{c}$, as defined below) at small $t$, $Q(t)\sim 1$ where as for large $t,Q(t)\sim 0$ signaling completion of structural relaxation. We evaluate $w(r)$ by monitoring the displacement of individual particles and once a particle moves by a distance greater than $r_{c}$ then we set $w(r)=0$ for all future times. This definition ensures $Q(t)$ is a monotonically decaying function of $t$.

Obviously, $Q(t)$ depends on the choice of $r_{c}$ and the optimal $r_{c}$ is chosen as follows: We study the fluctuation of $Q(t)$, defined as the dynamical susceptibility:

[TABLE]

at the lowest temperature for different choices of $r_{c}$. From definition, $\chi_{4}(t)=0$ when no particles move by a distance beyond $r_{c}$. It is zero, as well, when displacements of all particles are greater than $r_{c}$. So there must be a time scale, $\tau_{x}$, when $\chi_{4}(\tau_{x})$ reaches maximum for a given $r_{c}$. Thus, the maximum of $\chi_{4}$, i.e. the fluctuations in $Q(t)$, at $t=\tau_{x}$ signify the maximal dynamic heterogeneity. At a given $T$, we choose the value of $r_{c}$ for which such dynamic heterogeneity attains the maximum. For example, we find that for $N=150$, $\chi_{4}(t)$ reaches its maximum for $r_{c}=0.09$ which is about $14\%$ of $r_{0}$ (see Fig. 7(a)). At a given $T$, $\chi_{4}(t)$ attains maximum for a time proportional to the structural relaxation time Donati et al. (2002); Karmakar et al. (2014).

Fig. 7(b) shows the $t$-dependence of $Q(t)$ in irregular confinement at different $T$ for $N=150$. While at low $T,Q(t)$ decays very slowly, it falls to zero rapidly for $T>0.006$. The time dependence of $Q(t)\propto\mathrm{exp[-(t/\tau_{\alpha})^{c}]}$, beyond a very short time, is found exponential $(c=1)$ for $T>0.006$, whereas, it is stretched exponential $(c\approx 0.62)$ for $T=0.006$ and a much slower decay, seemingly power law for $T=0.002$.

To get an estimation of the structural relaxation time, $\tau_{\alpha}$, from the overlap function $Q(t)$, we use $Q(\tau_{\alpha})=1/{\mathrm{e}}$ for $T\geq 0.006$. For bulk system, one can estimate the structural relaxation time by looking into the temporal decay of the self part of the intermediate scattering function, $F_{s}(k,t)$ (with $k$ set to the value at which the static structure factor features its first peak). But, since the wave vector $k$ is not well defined for a finite system, we use $Q(t)$ to calculate $\tau_{\alpha}$.

Next, we study the $t$-dependence of $\chi_{4}(t)$ at different $T$ (Fig. 7(c)). We find that $\tau_{x}$ decreases to lower values with increase in $T$. We show the $T$ dependence of $\tau_{x}$ along with $\tau_{\alpha}$ in Fig. 7(d). The $T$-dependence of these two characteristic times appear to be identical within the tolerance. This implies that the enhanced heterogeneities cause structural relaxation in Coulomb clusters. The rapid increase of $\tau_{x}$ and $\tau_{\alpha}$ with decrease in $T$ are reminiscent of glassy systems.

IV.4.2 Cage-correlation time

The dynamics in Coulomb clusters can be further probed by addressing the ‘cage effects’. In glass formers, the particles get slower, without appreciable change in the structure, upon approaching the glass transition temperature. This is due to the caging effect in which each particle is locked up by its Coulomb repelling neighbours. Rearrangement of ‘cage’ relaxes the system and thereby the particles diffuse. The rate of change of surroundings of a particle yields a cage correlation function (CCF) and helps in understanding how rapidly the local environment of each particle changes on an average.

To calculate CCF we first define a neighbor list keeping track of each particle’s neighbors. If the list of a particle’s neighbors at time $t$ is identical to the list at time $t=0$, CCF assumes unity for that particle. If, however, any of the confined particle’s neighbor changes, CCF becomes zero at that instant of time. A neighbor list $\mathbf{L_{i}(t)}$ for particle $i$ in an $N$ particle system is a vector of length $N$, and is defined as Khademi et al. (2015)

[TABLE]

with $j=1,2,...,N$, and $f(r_{ij})=1$ if $j$ is the nearest neighbour of $i$ at time $t$ and zero otherwise. We use Voronoi construction to select neighbour list at all $t$. The CCF at $t$ is given by

[TABLE]

In Fig. 8(a) we show the decay of CCF at different $T$ for irregular confinement. We see that at $T=0.002$, $C_{g}(t)$ remains close to unity for all times, implying no significant rearrangements of neighbours. This is because, particles are just jiggling only around their equilibrium positions. As $T$ increases $C_{g}(t)$ shows a decay that turn out to stretched exponential in $t$: $c_{g}(t)\propto\mathrm{exp}\left[-t/\tau_{g}\right]^{\beta}$ with $\beta<1$ (details of the fitting parameters are given in Table I in supplementary materials). Such decay is found for both the confinements (for circular confinement see Fig. 7 in supplementary material). Similar decay is also common to supercooled liquids Rabani et al. (1999). Recently, it is found that water in nanoporous silica Shekhar et al. (2014) shows stretched exponential decay with the characteristics exponent $\beta=0.6$. Ref. Phillips (1996) describing stretched exponential temporal decay in supercooled liquids with $\beta\approx 0.6$. While the literature discusses only the short-range interacting systems, it is interesting to note that our Coulomb clusters (with long range interactions) show similar trends. For a system with with static random traps, the asymptotic value of $\beta$ relates the dimensionality $d$ of the system through $\beta=d/(d+2)$ Grassberger and Procaccia (1982). Thus for $d=2$, asymptotic value of $\beta$ is $0.5$. We find that for both the circular and irregular confinements $0.45<\beta<0.6$ which we consider is in broad agreement with such predictions.

We estimate the caging and non-caging times by studying the relative fluctuations in position of the particles with respect to its neighbors. A perfectly caged particle is expected to have equidistant neighbors. Exploiting this idea, we define the caging time $(\tau_{C})$ of the $i$-th particle as the time up to which the following condition holds for at least three of its neighboring particles $(j)$:

[TABLE]

Here $j=1,2,\cdots N_{b}$, $N_{b}$ denoting the number of nearest neighbors of particle $i$; $\Delta r_{ij}(t)=|\vec{r}_{i}(t)-\vec{r}_{j}(t)|$ and $\sigma^{2}_{i}(t)=\langle\Delta^{2}r_{ij}(t)\rangle-\langle\Delta r_{ij}(t)\rangle^{2}$. We take cut-off value $d$ as $0.1r_{0}$ (analogous to the usual Lindemann ratio). Once the above criterion breaks down, then the time taken by the particle to get caged again, is called the non-caging time $(\tau_{NC})$. In above, caging is defined at least with respect to three of its neighbors, because in $\mathrm{2D}$ a minimum of three particles are required to ‘cage’ a particle. Fig. 8(b) shows the averaged $\tau_{C}$ and $\tau_{NC}$ for irregular confinements. Here, $\tau_{C}$ and $\tau_{NC}$ are expressed relative to the total time of the MD simulation. At low $T$, caging time is expectedly high as particles are mostly confined by its nearest neighbors. As $T$ increases thermal energy overcomes caging. Thus, $\tau_{C}$ decreases with $T$ while $\tau_{NC}$ increases. From Fig.  8(b) we see that these two time scales crosses each other around the same $T_{X}$ that shows other signatures of the thermal crossover. We note that this analysis excludes the particles on the boundary because the current definition of ‘caging’ is ill defined for them.

IV.4.3 Persistence and exchange time

A crucial time scale that probes the relaxation mechanism as well as the glassiness of dynamics is the behavior of the persistence and exchange times Hedges et al. (2007); Berthier et al. (2005). They are defined as follows: At any given $T$, let us consider a particle $i$, whose initial $(t=0)$ position position is $r_{i}(0)$. The persistence time $t_{1}$ for a given cut-off distance $d$ is specified by the first time that particle $i$ has moved far enough that $|r_{i}(t_{1})-r_{i}(0)|\geq r_{p}$.

Exchange times $t_{n}-t_{n-1}$ for $n>1$, require the recursive determination $t_{n}$, setting $t_{n-1}$ as initial time. For a given $r_{p}$, distributions of exchange $(\tau_{e})$ and persistence $(\tau_{p})$ times can be obtained by ensemble averaging over all trajectories. $r_{p}$ is chosen of the order of $r_{0}$ so as to probe the structural relaxation and diffusion. The same $r_{p}$ is used for all $T$. It has been shown Hedges et al. (2007); Berthier et al. (2005) that for glassy dynamics, distribution of $(\tau_{e})$ and $(\tau_{p})$ generally decouple near the glass transition temperature. This decoupling is attributed to dynamic heterogeneity in these systems, and constitute a key feature of supercooled liquids Hedges et al. (2007); Berthier et al. (2005).

The distribution of these two time scales are shown in Fig. 9 for irregular (panel a) and circular (panel b) confinements with $N=500$. We notice that for $T>0.030$, $P(\tau_{e})$ and $P(\tau_{p})$ coincide for both the confinements. As $T$ becomes close to $T_{X}(\approx 0.020)$, the two distributions become distinct; mean value for $P(\tau_{p})$ moves toward longer times for IWM. In contrast, for circular confinement distribution of these two time scales do not appear to decouple. We have ensured that the decoupling is not because of bias in the sampling process by carrying out Jackknife test Press et al. (2007). Such decoupling is also observed for $N=150$ (see Fig. 8 in supplementary materials). Note that the distributions of $\tau_{e}$ and $\tau_{p}$ are shown for $T\geq T_{X}$ as at low $T$, displacement of particles beyond $r_{p}$ is rare and thus such distribution lacks confidence below $T_{X}$.

Since, the overlap function, cage correlation function and persistence time, all are connected with the structural relaxation of the system, we now look into the $T$-dependence of the characteristic time scales associated with these quantities.For example, $\tau_{g}$ (obtained by fitting the cage correlation function $c_{g}(t)$ as discussed) represents characteristic time scale for the rearrangement of particles in their local environment and average persistence time, $\langle\tau_{p}\rangle$, depicts the time required for particles to move a certain distance on average and thus associated with the relaxation mechanism in the system. In Fig. 10, we show the $T$-dependence of $\tau_{\alpha},\tau_{g}$ and $\langle\tau_{p}\rangle$ for irregular (panel a) and circular confinements (panel b). We find that all these time scales increase rapidly for IWM compared to CWM as $T$ approaches $T_{X}$. Thus, the signatures of glassy dynamics, while show up in both confinements, are more prominent for irregular confinement, broadly mimicking a bulk disorder system.

IV.4.4 Cross-Correlation between different relaxation times

Here, we study the cross-correlation between the different time scales described above. Fig. 10(c) illustrates the cross-correlation between $\tau_{\alpha}$ and $\tau_{g}$, $\tau_{\alpha}$ and $\langle\tau_{p}\rangle$, and $\langle\tau_{p}\rangle$ and $\tau_{g}$, gathered over a temperature window across melting; namely, $0.01\leq T\leq 0.05$, for irregular confinement. Our results indicate that there is a positive correlation between these characteristic times scale; when one increases, other also increases in a proportionate manner. Quantifying such correlations, we find that $\tau_{i}\propto\tau_{j}^{c_{ij}}$ with following sets of {$\tau_{i},\tau_{j},c_{ij}$}: {$\tau_{\alpha},\tau_{p},3.14\}$, $\{\tau_{\alpha},\tau_{g},3.98\}$ and $\{\tau_{p},\tau_{g},1.27\}$. Similar positive correlation between different times scales is also observed for circular confinement as depicted in Fig. 10(d). In this case, we find {$\tau_{i},\tau_{j},c_{ij}$}: $\{\tau_{\alpha},\tau_{p},0.73\}$, $\{\tau_{\alpha},\tau_{g},4.03\}$ and $\{\tau_{p},\tau_{g},4.45\}$. Thus, all these time scales give a coherent description of the nature of the dynamics of particles in traps.

V CONCLUSION

In conclusion, we have carried out comprehensive simulations to study the thermal crossover for Coulomb-particles in confinements, by probing their static and dynamic properties. While static properties are found quite insensitive to the symmetry of the confinements, dynamical correlations yield distinct signatures for the two confinements. We find intriguing stretched exponential spatial correlations in irregular confinements and a stretched Gaussian behavior in circular traps. Such a stretched Gaussian trend is recently observed in experiments with nano-particles in confined media He et al. (2014); Jacob et al. (2015). We devised several melting criteria, often taking advantage of the confined geometries, all of which point towards a single crossover temperature, when expressed in term of bulk-parameter $\Gamma$, yields $\Gamma\approx 147$ and and is quite close to the predicted critical value, $\Gamma_{c}\approx 137$, for the bulk system.

Our results offer insights into the ‘phases’ in a confined and finite system. From the study of static correlations, we find that the low-$T$ phase lacks positional order but posses bond-orientational order, which is further strengthened by temporal bond orientation correlation function showing an algebraic decay below the cross over temperature. We also witness key signatures of glassy behavior in our systems, such as, slow heterogeneous dynamics, rapid increase of relaxation time, decoupling of persistence and exchange times, etc. All these indicate that a hexatic-glass like phase is possible to realize in finite systems, at least in those with irregularity and perhaps for long-range inter-particle interactions. A hexatic glass, differing from the usual hexatic phase in $\mathrm{2D}$ only by its quenched defects – makes the study of defects in our system an interesting proposition. We hope that all these will help enriching our knowledge of thermal melting in two-dimension.

ACKNOWLEDGEMENTS

The authors acknowledge computational facilities at IISER Kolkata. BA acknowledges University Grant Commission (UGC), India, for doctoral fellowship. AG acknowledges the hospitality of the Aspen Center for Physics, and was supported in part by the Simons Foundation.

I Positional order in circular confinement

In this section, we illustrate and quantify the positional order in circular confinement with Coulomb interacting particles for a fair comparison with similar study in Fig. 1 of main article. The results in Fig. S.1(b) indicates that, in the ground state $(T=0)$, the positional order in $V_{\rm conf}^{\rm Cr}$ is better developed than in $V_{\rm conf}^{\rm Ir}$ (see Fig. 1 of main paper), producing more definitive peaks in the Fourier transform of density. This is even more prominent (Fig. S.1(c)) in the subsystem of selected particles where the stronger order is evident in naked eye (Fig. S.1(a)). At $T=0.020$, $\rho(\vec{k})$ shows a diffuse pattern (Fig. S.1(d)) similar to that observed for $V_{\rm conf}^{\rm Ir}$ (see Fig. 1(d)).

II Bond Orientational Order for circular confinement

Fig. S.2 shows the $T$-dependence of the distribution $P(|\Psi_{6}|)$ of bond orientational order parameter, $\left|\psi_{6}\right|$, for $V_{\rm conf}^{\rm Cr}$ with $N=500$ (panel a), $150$ (panel c) and $75$ (panel e) particles. $P(|\Psi_{6}|)$ shows a sharp peak around unity, depicting an orientationally ordered phase at low $T$, and the distribution becomes broader with increase in $T$. In Fig. S.2(b,e and f) we show the distribution $P(\phi_{6})$ of $\phi_{6}$ (see main text for definition), at different $T$ for $V_{\rm conf}^{\rm Cr}$ with $N=500,150$ and $75$, respectively. Here, we find that $P(\phi_{6})$ distinguishes the ordered solid-like phase at low-$T$ and orientationally disordered liquid-like phase at high-$T$. Thus, the temperature dependence of $P(|\Psi_{6}|)$ and $P(\phi_{6})$ for $V_{\rm conf}^{\rm Cr}$ are qualitatively similar to that discussed for $V_{\rm conf}^{\rm Ir}$ (see Fig. 2) in main paper. However, details vary with respect to $P(\phi_{6})$ for smaller systems (Fig. S.2(d,f)). $V_{\rm conf}^{\rm Cr}$ seems to present a gradual crossover for $N\leq 150$.

III Analysis of trajectory of particles in circular confinement

Fig. S.3 shows the trajectory of particles in circular confinement while corresponding analysis for irregular confinement is shown in Fig. 4 of the main paper. Motional features are qualitatively similar for particles in both the traps, but here, we emphasise following two points:

(i) Analysis of the instantaneous position $r_{i}(t)$ (Fig. S.3(d)) in $V_{\rm conf}^{\rm Cr}$ reveals that particles change their positions more gradually compared to those in $V_{\rm conf}^{\rm Ir}$ (see Fig. 4).

(ii) Circular symmetry dictates the development of coherent motion in $V_{\rm conf}^{\rm Cr}$. Lack of such spatial symmetry in case of $V_{\rm conf}^{\rm Ir}$ causes formation of tortuous string-like path through out the system (see Fig. 4).

IV Time evolution of exponent $k_{\rm Cr}$

In Fig. S.4 we show the time evolution of the exponent $k_{\rm Cr}$ (see main text for details), at different $T$ for $V_{\rm conf}^{\rm Cr}$ with (a)$N=500$, (b) $N=150$ and (c) $N=75$ particles. While at low $T(\approx 0.006)$ $k_{\rm Cr}$ continuously decays from $2$ to nearly $1$ (inset (a1), (b1) and (c1)) for all system sizes, $1\leq k_{\rm Cr}\leq 2$ for all higher $T$. Thus, the distribution of the displacement of the particles is stretched Gaussian for $V_{\rm conf}^{\rm Cr}$ while it is stretched exponential $(k_{\rm Ir}<1)$ for $V_{\rm conf}^{\rm Ir}$ (see Fig. 5) as discussed in the main article. Interestingly, stretched Gaussian spatial decay has been realized in recent experiments with nano-particles (see Ref. 80-81 in the main paper). Note that for $N=75$, $1\leq k_{\rm Cr}\leq 2$ but $k_{\rm Ir}<1$ (see Fig. 5(c)) for $T\geq 0.02$. Thus, standard liquid like behavior is not recovered for $V_{\rm conf}^{\rm Ir}$ with small number of particles.

V Temporal bond orientation correlation function for Circular confinement

While Fig. 6 in the main article shows the $t$-dependence of temporal bond orientational correlation function $(g_{6}(t))$ at different $T$ for $N=500$ particles, here in Fig. S.5 we discuss the time evolution of $g_{6}(t)$ for irregular (panel a) and circular (panel b) confinements with $N=150$ particles, respectively. For $N=150$, the $t$-dependence of $g_{6}(t)$ at different $T$ for both the traps is qualitatively similar to what is discussed in the main article. Here, we find that the algebraic decay of $g_{6}(t)$ persists over wider temperature window for $N=150$ (Fig. S.5(a)) compared to $N=500$ particles in $V_{\rm conf}^{\rm Ir}$ (see Fig. 6(a)). At intermediate $T$, data show larger fluctuations at long time for both the confinements due to smaller system size.

VI Overlap function $Q(t)$ for circular confinement

Fig. S.6(a) shows the $t$-dependence of the overlap function $Q(t)$ at different $T$ for $N=150$ particles in $V_{\rm conf}^{\rm Cr}$. At low $T$, the decay of $Q(t)$ is quite faster than that in irregular confinement (Fig. 7(b)) and shows stretched exponential decay. At high $T$, Q(t) decays exponentially to zero. Fig. S.6(b) shows $t-$dependence of $\chi_{4}(t)$ at different $T$. Here, we see that with increase in $T$, the time, $\tau_{x}$, at which $\chi_{4}(t)$ attains the maximum value decreases gradually (inset of Fig. S.6(b)). As discussed in the main paper, we have extracted the structural relaxation time $\tau_{\alpha}$ from the condition $Q(\tau_{\alpha})=1/e$. The $T$-dependence of $\tau_{\alpha}$ and $\tau_{x}$ are shown in the inset of Fig. S.6(b). While both the time scales show similar $T$ dependence, but for $V_{\rm conf}^{\rm Ir}$ these two time scales increase more rapidly with decrease in $T$ (see Fig. 7(d)). Thus, $T$-dependence of structural relaxation time indicates ‘better glassiness’ in $V_{\rm conf}^{\rm Ir}$.

VII Details of the fitting parameters for Cage correlation function

The cage correlation function, $c_{g}(t)$, is fitted with the following functional form:

[TABLE]

For both the confinements, $c_{g}(t)$ shows stretched exponential decay $(\beta<1)$ and values of $\tau_{g}$ and $\beta$ at different $T$ is given in Table. 1. Here, we find that $\tau_{g}$, the characteristic time scale associated with the local rearrangement of particles, grows more rapidly with decrease in $T$, for irregular confinement, reinstating the glassy dynamics in $V_{\rm conf}^{\rm Ir}$.

VIII Cage correlation function, caging and non-caging time

In Fig. S.7(a) we show the decay of cage correlation function with time at different $T$ for circular confinement. While at low $T$, $C_{g}(t)$ is close to unity, it shows stretched exponential decay at higher $T$. The details of the characteristic time associated with this decay is given in Table. 1. Fig. S.7(b) shows the $T$-dependence of average caging $(\tau_{C})$ and non-caging time $(\tau_{NC})$ for circular confinement. For $V_{\rm conf}^{\rm Cr}$, we find that the fraction of solid like particles is quite higher compared to $V_{\rm conf}^{\rm Ir}$ at low $T$ (see Fig. 8(b)). We also find that these two time scales cross each other around the crossover temperature, $T_{X}\approx 0.020$, as discussed in main paper.

IX Distribution of persistence and exchange times

Fig. S.8 shows the distribution of persistence and exchange times for irregular (panel a) and circular (panel b) confinements with $N=150$ particles. For $T>0.030$, the exchange $(\tau_{e})$ and persistence $(\tau_{p})$ time distributions coincide for both the confinements. As $T$ becomes close to $T_{X}(\approx 0.020)$, the two distributions become distinct for irregular confinement; mean value for the distribution of persistence time moves toward longer times (Fig. S.8(a)). Thus, we see the signature of decoupling, similar to what is described in main text for $N=500$, even for smaller systems. For circular confinement, while there is no sign of decoupling between $\tau_{e}$ and $\tau_{p}$, but for $T\geq 0.10$, all the distributions coincide as shown in Fig. S.8(b).

Thus, from our analysis we find that, the static properties of Coulomb interacting particles in circular and irregular confinements show qualitatively similar temperature dependence but the dynamical responses make clear distinction between the two. Signatures of slow and heterogeneous dynamics, while present in two systems, the qualitative measures differ substantially. Azimuthal symmetry of the $V_{\rm conf}^{\rm Cr}$ leads to circularly symmetric spatially heterogeneous dynamics of the particles at low temperatures.