---
abstract: 'We consider a two dimensional colloidal dispersion of soft-core particles driven by a one dimensional stochastic flashing ratchet that induces a time averaged directed particle current through the system. It undergoes a non-equilibrium melting transition as the directed current approaches a maximum associated with a resonance of the ratcheting frequency with the relaxation frequency of the system. We use extensive molecular dynamics simulations to present a detailed phase diagram in the ratcheting rate- mean density plane. With the help of numerically calculated structure factor, solid and hexatic order parameters, and pair correlation functions, we show that the non-equilibrium melting is a continuous transition from a quasi-long ranged ordered solid to a hexatic phase. The transition is mediated by the unbinding of dislocations, and formation of compact and string-like defect clusters.'
author:
- Shubhendu Shekhar Khali
- Dipanjan Chakraborty
- Debasish Chaudhuri
title: 'Structure-dynamics relationship in ratcheted colloids: Resonance melting, dislocations, and defect clusters'
---

§ ł Ł ø

\#1

Introduction {#sec:introduction}
============

A class of non-equilibrium driven systems called pump models are particularly intriguing due to their following property. They involve periodic forces, in time and space, that vanish under spatio-temporal averaging but still drives an overall directed current [@Julicher1997a; @Astumian2002; @Hanggi2009; @Reimann2002; @Brouwer1998; @Citro2003; @Jain2007; @Chaudhuri2011; @Chaudhuri2015; @Chaudhuri2015f]. This is achieved via the breaking of time-reversal symmetry through, e.g., a phase lag between spatially non-local drives [@Brouwer1998; @Jain2007; @Chaudhuri2011], or breaking of space inversion symmetry of the external potential profile [@Julicher1997a; @Reimann2002; @Astumian2002; @Hanggi2009]. Most of the biological processes generating directed motion involve reaction cycles and utilize some variant of this principle. Natural examples involve ion-pumps, e.g., the Na$^+$, K$^+$-ATPase pumps, and molecular motors [@Gadsby2009], e.g., Kinesin or myosin moving on polymeric tracks of microtubules or F-actins, respectively [@Reimann2002]. The flashing ratchet model has been used to describe molecular motor locomotion [@Julicher1997a]. In experiments on colloids, ratcheting could be generated using optical [@Faucheux1995; @Lopez2008], magnetic [@Tierno2010; @Tierno2012] or electrical fields [@Rousselet1994; @Leibler1994; @Marquet2002]. Most of the studies on pump models focused on systems of non-interacting particles, restricted to one dimension, with a few exceptions that analyzed the impact of interaction on molecular motors [@Derenyi1995; @Derenyi1996], collective properties of particle pumps [@Jain2007; @Marathe2008; @Chaudhuri2011; @Chaudhuri2015; @Chaudhuri2015f], and in ratchet models [@Savelev2004; @Pototsky2010; @Savelev2003; @Hanggi2009].

In a recent study, we used an asymmetric periodic potential that switches between an [*on*]{} and [*off*]{} state in a stochastic manner to drive a directed current of particles in a two dimensional (2d) dispersion of sterically stabilized colloids [@Chakraborty2014], focusing on the frequency and density dependence of the ratcheted current. With the change in the rate of ratcheting, the time- averaged directed current carried by the colloids show a resonance with the system’s relaxation frequency [@Chakraborty2014]. The current shows a non-monotonic dependence on density as well. This change in the dynamical properties, as we show in this paper, is closely related to the associated structural changes, e.g., the solid melts near the resonance frequency.

In the limit of extremely high switching frequency, higher than the inherent relaxation time of the colloids, the system can only respond to essentially a time- averaged potential profile. In addition, if one considers the limit of vanishing asymmetry in the potential profile, the scenario becomes equivalent to that of the re-entrant laser induced melting transition (RLIM) [@Chowdhury1985; @Wei1998; @Frey1999; @Chaudhuri2006], in which a high- density colloidal liquid undergoes solidification followed by melting, as the strength of a commensurate external periodic potential is increased. This is an equilibrium phase transition of the Kosterlitz-Thouless type [@Frey1999; @Chaudhuri2006], and is described in terms of unbinding of a specific type of dislocations, allowed by the potential profile.

In this paper we consider an asymmetric ratcheting of soft-core particles, and investigate structural transitions associated with the change in dynamical behavior of the system, observed in terms of its current carrying capacity. Using a large scale molecular dynamics simulation, we obtain the phase diagram in the density- ratcheting rate plane, showing melting from a solid to hexatic phase. We find a re-entrant solid- hexatic- solid transition with changing ratcheting frequency. The transitions are associated with a non-monotonic variation of the mean directed current. As we demonstrate in detail, the non-equilibrium melting is a continuous transition from a quasi- long ranged ordered (QLRO) solid to a hexatic phase, and is mediated by the formation of topological defects. The dominant defect types generated at the solid melting are dislocations, and compact or string-like defect clusters.

![image](phase_diag5.png){width="\linewidth"}

In \[sec:model\_simulation\] we present the model and details of numerical simulations. In \[sec:results\_discussion\] we discuss the detailed phase diagram, explaining the properties of the different non-equilibrium phases. The associated variation of driven directed current with driving frequency and density is shown in \[ssec:current\]. In this section we establish the relation of changing particle current to the non-equilibrium phase transitions. This is followed by a detailed analysis of the melting transitions in terms of the order parameters presented in \[ssec:reentrant\_transition\]. In the following three subsections, the phase- transitions are further characterized in terms of the distribution functions of order parameters, correlation functions, and formation of topological defects. We finally conclude presenting a discussion and outlook in \[sec:outlook\].

Model and Simulation Details {#sec:model_simulation}
============================

We consider a two dimensional system of a repulsively interacting colloidal suspension of $N$ particles in a volume $A=L_x L_y$. The mean inter-particle separation in this system $a^2=\sqrt{3} \rho/2$ is set by the particle density $\rho = N/A$. We assume that the colloids repel each other via a shifted soft-core potential $U(r)=\e \,[(\sigma/r)^{12}-2^{-12}]$ when the inter-particle separation $r< r_c r_c=2 \sigma$, and $U(r)=0$ otherwise. The units of energy and length scales are set by $\e$, $\s$ respectively. The system evolves under an asymmetric ratchet potential $U_{\rm ext}(x,y,t)=V(t) \left[\sin \left(2 \pi y/\lambda \right)+\alpha \sin \left(4 \pi y/\lambda \right) \right]$, where the time-dependent strength $V(t)$ switches between $\e$ and $0$ stochastically with a rate $f$. The two sinusoidal terms in the above expression of $U_{\rm ext}$ with $\a=0.2$ maintains the asymmetric shape of the potential profile. When it assumes a triangular lattice structure, the separation between consecutive lattice planes in the system is $a_y=\sqrt{3} a/2$. We have chosen the periodicity of the external potential $\lambda = a_y$, commensurate with the mean lattice spacing. In the absence of the external potential, the soft core solid is expected to undergo a two stage solid- hexatic- liquid transition [@Kosterlitz1973; @Halperin1978; @Young1979; @Kapfer:2015ca], with the solid melting point at $\r \s^2 \approx 1.01$. In the presence of a time- independent potential profile with $V(t)=U_0$ and $\a=0$, the system undergoes RLIM with increase in $U_0$ [@Wei1998; @Frey1999; @Chaudhuri2006]. At $U_0=\e$, the laser induced melting point of the soft-core solid is $\r \s^2 = 0.95$ [@Chaudhuri2006].

We perform molecular dynamics simulations of the system in the presence of an external ratcheting potential using the standard leap-frog algorithm [@Frenkel2002] with a time-step $\d t = 0.001\,\t$ where $\t= \s \sqrt{m/\e}$ is the characteristic time scale. We use $m=1$. The temperature of the system is kept constant at $T =1.0 \e/\kb$ using a Langevin thermostat characterized by an isotropic friction $\g = 1/\t$. At each step a trial move is performed to switch $V(t)$ between $0$ and $\e$, and accepted with a probability $f \, \d t$. In this paper, we present the results for a large system of $N=262144$ particles. We discard simulations over initial $10^7$ steps to ensure achievement of steady state, and the analyses are performed collecting data over further $10^7$ steps.

Results and Discussion {#sec:results_discussion}
======================

Phase diagram {#ssec:phase_diagram}
-------------

In \[fig:phase\_diagram\] we present the detailed phase diagram along with local density profiles of the 2D system of mono-dispersed ratcheted colloids. The system displays a solid and a density modulated hexatic phase, controlled by the dimensionless density $\r \s^2$ and ratcheting rate $f \t$. The color codes in \[fig:phase\_diagram\]($a$) denote the values of mean solid order parameter. The two dashed lines with open symbols signify the two solid melting boundaries at small and high frequencies. The dash-dotted line with filled squares denotes the inverse of relaxation time-scales $f_s$ of the system at a given density. For ratcheting rates slower than this time-scale, the solid and hexatic order can [*equilibrate*]{} to the instantaneous external potential and follow its change. The details of the calculation of such relaxation times are discussed in \[appendix:relaxtion\]. We characterize the various phases using the structure factor, the pair correlation function, the solid and hexatic order parameters, and their distribution functions. Further details of such characterization are presented in later sections.

As is shown in \[fig:phase\_diagram\]($a$), the system remains in a QLRO triangular lattice solid phase at the highest frequencies, if the ambient density permits. As the frequency decreases, the solid melts into a hexatic phase, below the dashed line through open circles (\[fig:phase\_diagram\]). As is shown later, the melting is a continuous transition and happens via a proliferation of topological defects, including the unbinding of dislocation pairs. The molten phase is a hexatic displaying a unimodal distribution of local- hexatic order parameter. This excludes any possibility of phase coexistence, indicating a continuous transition. As the frequency of the drive is decreased further, below the equilibrium relaxation time, the system starts to follow the time variation of the external potential. As a result, the time-averaged properties turn out to be approximately a superposition of the properties of the equilibrium states in the presence ($V(t)=\e$) and absence ($V(t)=0$) of external potential profile. At high densities ($\r \s^2 \gtrsim 1.026$) the solid phase is stabilized even at low ratcheting frequencies like $f\t=0.01$. The melting boundary of this solid is shown by the open $\triangledown$ and dashed line.

In \[fig:phase\_diagram\]($b$)-($h$) we present the time-averaged local density profiles $\r(x,y)$ at two mean densities $\r \s^2=1.04,\, 0.98$. We show the results over local cross- sections of area $10\s \times 10\s$, for better visibility. At $\r \s^2=1.04$, the density profile shows triangular lattice structure at $f \t=0.01$ (\[fig:phase\_diagram\]($b$)). At $f \t=1$ the solid melts into a phase with density modulation along $y$-axis, the direction of ratcheting drive (\[fig:phase\_diagram\]($c$)). The melting is quantified later in terms of vanishingly small solid order. At $f \t=10$, a local triangular lattice- like pattern reappears, albeit the solid order remains small (\[fig:phase\_diagram\]($d$)). The system shows a more compact triangular lattice solid at $f\t=30$ associated with appearance of significant solid order (\[fig:phase\_diagram\]($e$)). As is turns out, this order decreases with the system size in a power law manner, identifying a QLRO solid (see \[fig:scaling\_solid\_order\_parameter\]). As we show later, the density modulated phases of the system at $f \t=1$ and $10$ share similar amount of hexatic order (see \[fig:solid\_op\_freq\_density\]).

In the top and left axes around the $\r(x,y)$ plots, we show the linear density profiles $\r_x$ (ordinate) and $\r_y$ (abscissa) measured along the horizontal and vertical dashed white lines indicated in the $\r(x,y)$ plots. The clean density modulations in $\r_x$ captures the spontaneous emergence of the solid- like order in that direction. Due to the shape of the triangular lattice solid, $\r_y$ shows large followed by tiny peaks along the white line perpendicular to the lattice planes, e.g., in \[fig:phase\_diagram\]($b$) and ($e$). Note that this feature is not shared by the other local density plots after the solid melts, shown in \[fig:phase\_diagram\]. The strong density modulation $\r_y$ in the hexatic phase is induced directly by the external potential minima. This phase lacks density modulation in $\r_x$, as is shown in \[fig:phase\_diagram\]($c$). However, as we shall quantify later, it shows significant hexatic orientation order.

Similar characterization at the lower density $\r \s^2 = 0.98$ are presented in \[fig:phase\_diagram\]($f$)-($h$). A clear solid- like triangular lattice order is observed at $f\t=30$ (\[fig:phase\_diagram\]($h$)). \[fig:phase\_diagram\]($g$) shows local weak triangular lattice like pattern that gets smeared within the hexatic phase. As we show later in Sec.-\[ssec:reentrant\_transition\], this phase in \[fig:phase\_diagram\]($f$) and ($g$) indeed display a vanishingly small solid order, in the presence of a reasonably large hexatic order. The system remains in the molten phase at the intermediate frequency range, between $0.1 \leq f \tau \leq 10$. Note that the densities considered here are relatively large with respect to the equilibrium melting point in the absence of external potential. As we demonstrate in the following, the non-equilibrium melting discussed in this section is associated with an increase in directed particle current carried by the system under the flashing ratchet drive.

Direct particle current out of alternating drive {#ssec:current}
------------------------------------------------

The inherent asymmetry of the ratchet potential, added with the stochastic switching drives a time- and space- averaged directed particle current in the system, $$\label{eq:current_def}
\left\langle  j_y
\right\rangle=\frac{1}{\tau_m}\frac{1}{A}\int^{\tau_m}
\upd_t\int^{L_x}\upd_x\int^{L_y}\upd_y ~j_y(x,y,t),$$ where, $\tau_m$ is an integral multiple of $1/f$, the mean switching time of the flashing ratchet potential.

![(color online ) [(a)]{} Variation of the scaled particle current along the direction of the drive as a function $f/\nu$ for different densities as indicated in the legend. The solid line is plot of the function $(f/\nu)/\left(1+(f/\nu)^2\right)$ as given by \[eq:current\]. [(b)]{} Plot of the current amplitude $\rho \kappa v_0$ as a function of density. The solid line is a fit to the data using the functional form $\kappa D_0 \rho^{3/2}(1-\rho/\r_c)$ with $\kappa D_0$ and $\rho_c$ as fitting parameters with values $\kappa D_0 \approx 0.22$ and $\rho_c \s^2 \approx 1.004$. []{data-label="fig:fluxplot"}](flux_amp_combi.pdf){width="\linewidth"}

The competition between the intrinsic relaxation of the system and the external drive leads to a resonance in the mean particle current. At low frequencies the particle current increases linearly with the driving frequency $f$, achieves a maximum around $f=\nu$, and beyond this decays as $f^{-1}$. The behaviour of the current in the whole frequency range can be captured by the simple ansatz [@Chakraborty2014] $$\label{eq:current}
  \left\langle j_y \right\rangle = \kappa \frac{\nu f}{\nu^2+f^2}
  \rho v_0,$$ where, $\nu$ is the intrinsic relaxation frequency, $v_0 = \nu \l$ intrinsic velocity, and $\kappa$ is a proportionality constant.

In the high density regime, where, the mean-free path of the particles is small, the diffusive time scale $\tau_D$ to travel the typical distance $\lambda$ is set by $\tau_D=\lambda^2/D(\rho)$. Here $D(\rho)$ denotes the density- dependent tagged particle diffusivity, which we assume to decrease linearly with density, $D(\rho)=D_0(1-\rho/\rho_c)$ [@Lahtinen2001; @Falck2004]. Here $D_0=\kb T/\gamma$ is the bare diffusivity. The commensurate external potential ensures that $\lambda^2 \sim 1/\rho$, and consequently the intrinsic relaxation frequency takes the form $\nu=\rho D(\rho)$. The intrinsic velocity scale is set by $v_0=\lambda/\tau_D= \rho^{1/2}D(\rho)$. Substituting for $v_0$ and $\nu$ in \[eq:current\], the density and frequency dependent current takes the form $$\langle j_y \rangle =\frac{f D_0^2}{D_0^2 \rho^2 (1-\rho/\rho_c)^2+f^2} \rho^{5/2} (1-\rho/\rho_c)^2.
    \label{eq:current_expression}$$ The resonance in the particle current appears at the ratcheting rate $f=\nu=D_0 \rho (1-\rho/\rho_c)$ and the density dependent amplitude take the form $\kappa \rho v_0 = \kappa D_0 \rho^{3/2}(1-\rho/\rho_c)$.

![(color online ) Plot of the mean-square displacement at density $ \rho \s^2 = 1.04$ for frequencies $f \tau=0.01$ (green), $1.0$ (red) and $30.0$ (blue). The displacements in $\la \D x^2 \ra$ are shown with open symbols and that in $\la \D y^2 \ra$ are represented with filled symbols. Typical trajectories of particles at densities $\rho\sigma^2=1.04$ ($b$) and $0.98$ ($c$) are shown for three different frequencies $f\tau=0.01$, $1$ and $30$. At high frequencies the particles get localized, whereas at very low frequencies localization happens only at densities above $\rho_c$. In the intermediate frequencies, at both densities the trajectories span a length scale of $\approx 20\,\s - 30\,\sigma$. []{data-label="fig:msdplot"}](msd_rho104 "fig:"){width="0.7\linewidth"} ![(color online ) Plot of the mean-square displacement at density $ \rho \s^2 = 1.04$ for frequencies $f \tau=0.01$ (green), $1.0$ (red) and $30.0$ (blue). The displacements in $\la \D x^2 \ra$ are shown with open symbols and that in $\la \D y^2 \ra$ are represented with filled symbols. Typical trajectories of particles at densities $\rho\sigma^2=1.04$ ($b$) and $0.98$ ($c$) are shown for three different frequencies $f\tau=0.01$, $1$ and $30$. At high frequencies the particles get localized, whereas at very low frequencies localization happens only at densities above $\rho_c$. In the intermediate frequencies, at both densities the trajectories span a length scale of $\approx 20\,\s - 30\,\sigma$. []{data-label="fig:msdplot"}](traj_new "fig:"){width="\linewidth"}

\[fig:fluxplot\]($a$) shows data collapse of particle currents when plotted as $(\kappa \rho v_0)^{-1} \langle j_y \rangle $ against the dimensionless variable $f/\nu$. The current maximizes at the resonance frequency of $f=\nu$. \[fig:fluxplot\]($b$) shows the limit of validity of the approximate form $v_0 = D_0 \rho^{1/2}(1-\rho/\rho_c)$. A comparison of \[fig:phase\_diagram\] and \[fig:fluxplot\]$a$ shows the relationship between the structure and dynamics. For example, at the resonance frequency, the system melts in order to carry the largest directed current.

The dynamics at the local scale can be better appreciated by examining trajectories of individual particles. First we consider the system at high density $\rho \sigma^2=1.04$. Both the components of the displacement fluctuations in \[fig:msdplot\]($a$) show an initial ballistic part $\sim t^2$ due to the inertial nature of the dynamics. Both in the low and high frequency limits, this crosses over to a sub-diffusive regime ($\sim t^\nu$ with $\nu<1$). The test particles show localized motion within cages formed by neighbours (see \[fig:msdplot\]($b$)). The intermediate frequency driving at $f\t=1$ shows a long time diffusive behavior in $x$, $\la \D x^2 \ra \sim t$, and almost ballistic motion in $y$, $\la \D y^2 \ra \sim t^\nu$ with $1 < \nu \lesssim 2$ (\[fig:msdplot\]($a$)). The corresponding particle trajectories display long excursions in the $y$-direction, as is displayed in \[fig:msdplot\]($b$).

\[fig:msdplot\]($c$) shows typical particle trajectories at low densities, $\r \s^2 = 0.98$. They get localized only at the highest driving frequencies $f \t=30$. At small and intermediate frequencies they show long excursions which are extended mainly in the $y$-direction, the direction of ratchet drive. The localization of trajectories is related to the low directed current carried by the system, and its ordering into a solid phase. Similarly, the extended particle trajectories are related to the melting of the solid and the presence of relatively large directed current.

![Plots of the static structure factor $\la \psi_{\bf q} \ra$ for the densities $\rho \sigma^2=1.04$ ($a$-$c$) and $0.98$ (figures $d$-$f$). The three columns correspond to three different frequencies $f \tau=0.01$ ($a$ and $d$), $1.00$ ($b$ and $e$) and $30$ ($c$ and $f$).[]{data-label="fig:struc"}](skplane_driven_sys_plot.png){width="8.6cm"}

Non-equilibrium melting {#ssec:reentrant_transition}
-----------------------

For a more quantitative analysis, we turn our attention to the phase diagram \[fig:phase\_diagram\]($a$) and consider the phase behavior along constant frequency, and constant density lines. The structure factor, $\la \psi_{\bf q} \ra = (1/N)\,\la \rho_{\bf q}\, \rho_{-\bf q} \ra$ (see \[fig:struc\]) with $\rho_{\bf q}=\sum_{j=1}^{N} e^{-{\mathbbm{i}}{\bf q} \cdot \bf{r}_j}$ and $\rho^*_{\bf q}=\rho_{-\bf q}$, can clearly distinguish between a solid, hexatic, liquid and a modulated liquid phase [@Chaikin2012]. In the solid phase $\la \psi_{\bf q} \ra$ shows a characteristic six fold symmetry with peaks at $\G_1=(0,\pm 2 \pi /a_y)$ and $\G_2=(\pm 2\pi/a,\pm \pi/a_y)$ reflecting the underlying triangular lattice structure (see \[fig:struc\]($a$),($c$),($f$)). The six intensity maxima broaden along the constant radius $q=2\pi/a$ circle in a hexatic (\[fig:struc\]($b$),($e$)). In a simple liquid with spherical symmetry, the broadening extends to overlap forming a characteristic ring structure. On the other hand, in a modulated liquid phase, $\la \psi_{\bf q} \ra$ is expected to show two bright spots at $\G_1$, in addition to the ring structure characterizing a simple liquid (see \[fig:struc\]($b$),($e$)). Note that the presence of the external periodic potential in the present context induces an explicit symmetry breaking by imposing density modulations in $y$-direction, ensuring $\la \psi_{\G_1} \ra> \la \psi_{\G_2} \ra$. The other four quasi Bragg peaks at $\G_2$, e.g., at the highest ratcheting frequencies, identify the appearance of the quasi long ranged positional order (QLRO). We use their arithmatic mean as the measure of solid order parameter $\la \psi_{\G_2} \ra$.

The QLRO in the solid phase is explicitly demonstrated using the system size dependence of $\la \psi_{\G_2} \ra$ shown in \[fig:scaling\_solid\_order\_parameter\]. The calculations are performed over sub-blocks of sizes $\ell_x \times \ell_y = \ell^2$. Within both the high and low frequency solid regions, $\la \psi_{\G_2}(\ell) \ra \sim \ell^{-\nu}$ where $\nu < 1/3$, the value of the exponent expected at the equilibrium KTHNY (Kosterlitz- Thouless- Halperin- Nelson- Young) melting [@Kosterlitz1973; @Halperin1978; @Young1979]. In this case, the value of the exponent $\nu$ depends on the mean density and ratcheting frequency.

![System size dependence of the solid order parameter. The ratio of mean order parameter measured over blocks of size $\ell= (\ell_x \times \ell_y)^{1/2}$ with respect to that measured over the whole system, $\la \psi_{\mathbf{G}_2}(\ell)\ra/\la \psi_{\mathbf{G}_2}(L) \ra$ decays with $\ell/L$ with power law $(\ell/L)^{-\nu}$. The data is shown for different densities as indicated in the legend. The solid and open symbols denote results at frequency $f\t=30$ and $f\t=0.01$, respectively. The solid black line is a plot of the power-law $(\ell/L)^{-1/3}$ expected at the equilibrium KTHNY melting point.[]{data-label="fig:scaling_solid_order_parameter"}](psiG2_vs_boxL_combi_3.pdf){width="6cm"}

![ Dependence of the ($a$)solid order parameter $\langle \psi_{\bf{G}_2} \rangle $ and ($b$)hexatic order parameter $\langle \psi_6 \rangle$ as a function of frequency is plotted for different densities as indicated in the legend on the top. The equilibrium melting point $\la \psi^{m}_{\bf{G}_2} \ra = 0.376$ is denoted by the dashed lines in ($a$). The variations of ($c$)$\langle \psi_{\bf{G}_2} \rangle$ and ($d$)$\langle \psi_6 \rangle$ are shown as a function of the mean density for three driving frequencies $f\tau = 0.01$, $1$ and $30$ as indicated in the legends. []{data-label="fig:solid_op_freq_density"}](op_combi_2.pdf){width="\linewidth"}

The phase behaviors are further characterized by following the change in the hexatic bond orientational order $\la \psi_6 \ra = (1/N) \la | \sum_{i=1}^N \ \vec \psi^i_6|^2 \ra$, where we define the local hexatic order $\vec \psi^i_6 = \sum_{k=1}^{n_v} (\ell_k/\ell) \exp({\mathbbm{i}}6 \phi_{ik}
)$ utilizing the $n_v$ Voronoi neighbors of the $i$-th test particle. Here $\phi_{ik}$ is the angle subtended by the bond between the $i$-th particle and its $k$-th Voronoi neighbor. In this definition we used the weighted average over the weight factor $\ell_k/\ell$ such that $\ell = \sum_{k=1}^{n_v} \ell_k$ and $\ell_k$ denotes the length of the Voronoi edge corresponding to the $k$-th topological neighbor [@Mickel2013].

In \[fig:solid\_op\_freq\_density\]($a$) and ($b$) we show the variations of the solid and hexatic orders, $\langle \psi_{\G_2} \rangle$ and $\langle \psi_6 \rangle$, as a function of the ratcheting frequency keeping the density of the system fixed. \[fig:solid\_op\_freq\_density\]($c$) and ($d$) show similar plots, but as a function of the density, keeping the driving frequency fixed. Both the order parameters show non-monotonic variation with frequency. The large variation of the solid order parameter $\la \psi_{\G_2} \ra$ with $f\t$ signifies melting, followed by a [*re-entrant solidification*]{}. Here we use the value of the solid order parameter at the equilibrium melting point, $\la \psi^m_{\G_2} \ra \approx 0.376$ (\[appendix:PT\]), to identify the boundary of solid phase, $\la \psi_{\G_2}\ra \ge \la \psi^m_{\G_2} \ra$. As we show later, the density-density correlation changes from a power law to an exponential decay across this melting.

As \[fig:solid\_op\_freq\_density\]($a$) shows, at all the densities considered the system remains in a solid phase at the highest driving frequencies. With reduction of frequency below $f\t \approx 10$, the system melts. At densities $\r \s^2 \gtrsim 1.01$ the system re-solidifies as the frequency is lowered further below $f\t \approx 0.1$. At intermediate to high frequencies, whether the system remains in a solid or fluid phase is essentially determined by the driving frequency, and not the mean density of the system (see \[fig:solid\_op\_freq\_density\]($c$)). Only at the lowest driving frequencies, one finds density dependent solidification. The hexatic order parameter $\la \psi_6\ra$ shows similar variations, albeit with lesser magnitude (see \[fig:solid\_op\_freq\_density\]($b$), ($d$)). At driving frequencies much larger than the inverse relaxation times, the system responds to only a time-integrated potential profile. The corresponding behavior is similar to that in the presence of a time-independent commensurate potential. For equilibrium melting of the free system, and melting in the presence of a time-independent periodic potential commensurate with the density, see \[fig:op\_free\_lif\] presented in \[appendix:PT\].

![Probability distributions of the local solid and hexatic order parameters $\psi^l_{\bf{G}_2}$ and $|\vec \psi_6^i|^2$. The distribution functions ${\cal P}(\psi^l_{\bf{G}_2})$ at densities $\rho \sigma^2=0.98$($a$) and $\rho \sigma^2=1.04$($b$). The local solid order parameter is determined over subsystems of length $\ell = (\ell_x \times \ell_y)^{1/2}$, where $\ell_x/L_x=\ell_y/L_y=1/14$. The distribution functions of local hexatic order ${\cal P}(|\vec \psi_6^i|^2)$ at $\rho \sigma^2=0.98$($c$) and $1.04$($d$), corresponding to four representative ratcheting frequencies as indicated in the legends.[]{data-label="fig:solid_op_hist_sq"}](lpsiG2_hist_driven_sys_combi.pdf "fig:"){width="\linewidth"} ![Probability distributions of the local solid and hexatic order parameters $\psi^l_{\bf{G}_2}$ and $|\vec \psi_6^i|^2$. The distribution functions ${\cal P}(\psi^l_{\bf{G}_2})$ at densities $\rho \sigma^2=0.98$($a$) and $\rho \sigma^2=1.04$($b$). The local solid order parameter is determined over subsystems of length $\ell = (\ell_x \times \ell_y)^{1/2}$, where $\ell_x/L_x=\ell_y/L_y=1/14$. The distribution functions of local hexatic order ${\cal P}(|\vec \psi_6^i|^2)$ at $\rho \sigma^2=0.98$($c$) and $1.04$($d$), corresponding to four representative ratcheting frequencies as indicated in the legends.[]{data-label="fig:solid_op_hist_sq"}](lpsi6_hist_all_snaps_combi.pdf "fig:"){width="\linewidth"}

![ (Color Online) The pair correlation functions $g(x,0)-g(\infty,0)$ calculated at the mean densities $\rho \sigma^2=0.98$ ($a$) and $1.04$ ($b$) for the frequencies indicated in the legends. []{data-label="fig:gx_plot"}]({gx_combi_3}.pdf){width="0.65\linewidth"}

Continuous transition: Distribution of order parameters
-------------------------------------------------------

We probe the order of the phase transitions using the distribution of the local solid and hexatic order parameters. In determining the local solid order, we divided the simulation box into sub-boxes of size $\ell = (\ell_x \times \ell_y)^{1/2}$ with $\ell_x/L_x=\ell_y/L_y=1/14$. $\psi^l_{\G_2}$ is then calculated using the definition of $\psi_{\bf q}$ restricted within these sub-boxes. For the local hexatic order parameter, we calculate $|\rho^i_6|^2$ for all the particles. The distribution functions of these quantities, ${\cal P}(\psi^l_{\G_2})$ and ${\cal P}(|\rho_{6}|^2)$ are plotted in \[fig:solid\_op\_hist\_sq\]. They remain unimodal at all points of the phase diagram. At low densities, $\r \s^2 = 0.98$ in \[fig:solid\_op\_hist\_sq\]($a$), the maximum of ${\cal P}(\psi^l_{\G_2})$ appears at an order parameter corresponding to the solid phase only at the highest frequencies. With decreasing frequency the peak shifts towards lower values, signifying melting below $f\t = 10$ and remain low at the lowest frequencies. The unimodal nature of the distribution function, as the peak shifts to lower values, signifies the absence of any metastable state across the transition, a characteristic of continuous transitions. At high densities, e.g., $\r \s^2 = 1.04$ in \[fig:solid\_op\_hist\_sq\]($b$), corresponding to the re-entrant transition, the peak of the distribution ${\cal P}(\psi^l_{\G_2})$ shifts from a high to low and back to high values as the ratcheting rate decreases from the highest frequencies. As before, the unimodal nature of the distribution corresponds to a continuous melting transition.

The solid melts to a hexatic phase, characterized by the finite hexatic order. However, a further melting of the hexatic is not observed as the frequency is varied. In the density-frequency range bounded between the two dashed lines with open inverted triangles and open circles denoted in \[fig:phase\_diagram\] ($a$), the system remains in a hexatic phase. This is corroborated by the distribution of the local hexatic order ${\cal P}( |\rho_6|^2)$ shown in \[fig:solid\_op\_hist\_sq\]($c$) and ($d$) corresponding to densities $\rho \sigma^{2}=0.98$ and $1.04$, respectively. The uni-modal nature of the distribution with a roughly unchanged peak position and a fat tail persists throughout the frequency range. The peak of the distribution does not shift. However, it is important to note that deep inside the hexatic phase, near $f \t = 1$, a significant fraction of the system displays vanishing hexatic order. This is more prominent at lower densities (see \[fig:solid\_op\_hist\_sq\]($c$)). As we show in a later section, local dip in the hexatic order is associated with the formation of grain boundaries.

Melting of solid: Correlation functions {#ssec:pair_correlation}
---------------------------------------

The pair correlation functions $g(x,y)$ capture the solid melting (see \[appendix:pair\_correlation\_function\] for further details). A density modulation is externally induced in the system along the $y$-direction by the ratcheting drive, breaking the translational symmetry in that direction, explicitly. To study the spontaneous symmetry breaking, here we focus on the $x$-component of the two point correlation functions $g(x,y)$. The component of the correlation $g(x,0)-g(\infty,0)$ along the axis perpendicular to the direction of the ratcheting drive are shown in \[fig:gx\_plot\]. This provides a more conclusive evidence to the nature of the phases in different density and frequency regimes. At low densities, as we have shown before, the solid order exists only at the highest frequencies. At such frequencies we find an algebraic decay of the correlation along the $x$ direction, signifying the QLRO (see \[fig:gx\_plot\]($a$)). At intermediate and low frequencies the system melts. This is captured by the exponential decay of the correlation with correlation length $\sim 10\sigma$ (see \[fig:gx\_plot\]($a$)). The scenario changes at higher densities. At high frequency, as before, we again find a solid phase, with the correlation exhibiting algebraic decay corresponding to the QLRO (\[fig:gx\_plot\]($b$)). At high densities, one obtains another solid phase at the low ratcheting frequencies. This also shows algebraic decay of correlations signifying QLRO (see $f\t =0.01$ graph in \[fig:gx\_plot\]($b$)). The power law $x^{-1/3}$ shown in the figures denote the expected correlation at the onset of the equilibrium KTHNY melting. At the intermediate driving frequencies, $f\t=1$ in $\r\s^2=1.04$ system, the correlation shows exponential decay with a correlation length $\sim 10\,\s$, similar to the behavior observed in the low density regime (\[fig:gx\_plot\]($b$)). The change from algebraic to exponential decay is utilized to identify the solid melting points in the phase diagram, as is detailed further in \[appendix:pair\_correlation\_function\].

Defect formation
----------------

The equilibrium melting of the two- dimensional solid within the continuous KTHNY melting scenario is known to proceed by the unbinding of dislocation pairs into free dislocations. To identify such topological defects, we first obtain the coordination number $n_v$ of each particle in the system counting the number of its Voronoi neighbours. In a perfect triangular lattice $n_v=6$ for all particles. We follow $n_v \neq 6$ particles to identify the $n_v$-fold defects. Even within the solid phase, fluctuations of bound quartets of $5-7-5-7$ defects (bound dislocation pairs) keep appearing. They form dislocations by dissociating into separate $5-7$ and $7-5$ non-neighboring pairs. Presence of a finite fraction of particles associated with dislocations characterize the hexatic phase. The system shows dislocation formation as the solid melts. Moreover, we find defect clusters larger than quartets that are either compact or string-like (grain boundary) [@Qi:2014]. All the dominant defect types observed in our simulations are indicated in \[fig:defect\_types\]($f$). Their typical configurations in a sub-volume of size $100 \s \times 100 \s$ at $\r \s^2=1.04$ and different ratcheting frequencies are shown in \[fig:defect\_types\]($a$)-($d$). In these figures, the colors associated to particles indicate the number of topological neighbors they have, $n_v=4\,$(purple), $5\,$(red), $6\,$(green), $7\,$(blue), $8\,$(yellow). Clearly, defect formation is suppressed at both the extremities of the ratcheting frequency. It increases significantly in the intermediate frequency regime associated with solid melting (\[fig:defect\_types\] ($b$)). The relative fraction of different defect types also vary with the driving frequency. In the highest frequency solids, only bound quartets (bound dislocation pairs) are observed in \[fig:defect\_types\] $(d)$. As the solid melts with decreasing frequency, dislocations and defect clusters start to appear and eventually dominate over the quartets in the system (\[fig:defect\_types\] $(b)$ and ($c$)). The string-like defects remain extended along the $y$-direction, the direction of particle current under ratcheting. Such a connected string of defects is shown in \[fig:defect\_types\] $(b)$ and has been highlighted in \[fig:defect\_types\] $(e)$, which shows Voronoi diagram of a region containing the connected string of defects. The color code in each Voronoi cell denotes the amount of hexatic order, and the arrows denote the corresponding hexatic orientations. At the location of the connected clusters of defects, the local hexatic order is low, and shows a hexatic orientation approximately orthogonal to the neighboring defect-free regions. With further lowering of the ratcheting frequency below $f_s$, the defect fraction decreases strongly.

This description is quantified by focusing on the time evolution of defect fractions. We first consider the evolution of the total fraction of all the topological defects, the percentage of particles having non-six Voronoi neighbors $N_d = (1 - n_6/N)\times 100$, where $n_6$ denotes the total number of particles with $n_v=6$ (\[fig:defect\_dynamics\] $(a)$). Clearly, the largest value of $N_d$ with the strongest fluctuations appear at the intermediate frequencies. The defect formation gets dramatically suppressed in the solid phase corresponding to the high ratcheting frequencies. At the lowest frequencies ($f\t=0.01$ in \[fig:defect\_dynamics\] $(a)$), $N_d$ remains relatively low and follows the switching of external potential.

The mean value $\la N_d \ra$ remains less than $4\%$ and varies non-monotonically with $f\t$ (\[fig:defect\_dynamics\] $(b)$). It shows a maximum at the resonance frequency corresponding to the largest directed current, relating formation of topological defects with carrying capacity of particle current in the system.

Further insight into the structure- dynamics relations can be obtained by following the behavior of different defect fractions separately. For this purpose, the percentage fraction $d_f$ of a defect type is defined as $d_f=(n_d/N) \times 100$, where $n_d$ is the total number of $n_v\neq 6$ particles that may contribute to either a quartet, a dislocation, a cluster, or a disclination as described above. The time averaged percentage fractions of these topological defects as a function of the driving frequency is shown in \[fig:defect\_dynamics\] $(c)$. They exhibit a similar non-monotonic behavior as $N_d$ and the mean particle current. In the high frequency solid, the dominant defects are the quartets. As the frequency is decreased, the melting of the solid is mediated by the unbinding of these quartets into dislocations. The dislocation fraction becomes larger than that of quartets. More importantly, at the resonance melting, the formation of defect clusters dominate (\[fig:defect\_dynamics\]$(c)$). The fraction of disclinations remain relatively insignificant (less than $0.03\%$), about two orders of magnitude smaller than that of the defect clusters. This is consistent with the fact that the hexatic does not melt within these parameter regimes.

Discussion {#sec:outlook}
==========

In conclusion, using a large scale simulation involving $262144$ particles, we have presented a detailed study of a ratcheted two-dimensional colloidal suspension, focusing on the structure- dynamics relationship. The mean directed particle current driven by the ratchet exhibits a resonance behavior. Associated with this, the solid melts to hexatic providing a mechanism allowing directed transport. The system exhibits a rich non-equilibrium phase diagram as a function of the driving frequency and mean density. At high densities, we found a re-entrant melting transition as a function of ratcheting frequency. The different phases are characterized by the spatially resolved density profile, the density-density correlation function, the structure factor, the solid and hexatic order parameters, and their distribution functions. The role of the defects in the phase transition has been investigated in detail. The solid- melting is associated with formation of dislocations, but unlike the equilibrium two- dimensional melting, this non-equilibrium melting is dominated by the formation of defect clusters, connected strings of defects that remain oriented largely along the direction of the ratcheting drive. Remarkably, the driven hexatic does not melt to a fluid within the studied range of density and ratcheting drive. Our detailed predictions regarding the variation of particle current and associated phase transitions can be verified using colloidal particles and optical [@Faucheux1995] or magnetic ratcheting [@Tierno2012] in a suitable laser trapping setup [@Wei1998]. The impact of changing degree of potential asymmetry on the dynamics and phase behavior remains an interesting future direction of study.

Acknowledgements {#acknowledgements .unnumbered}
================

DC thanks ICTS-TIFR, Bangalore for an associateship, and SERB, India, for financial support through grant number EMR/2016/001454.

Relaxation under external potential {#appendix:relaxtion}
===================================

![Relaxation dynamics of the solid and hexatic order parameters at $\r \s^2 =1.04$($a$) and $\r \s^2 =0.98$($b$), respectively. The time scales for relaxation are $f_s^{-1} \approx 6\,\t$($a$) and $f_s^{-1}\approx 70\,\t$($b$).[]{data-label="fig:relaxation_frequency"}](relaxation_combi_new2.pdf){width="50.00000%"}

![Plot of the solid order parameter $\la \psi_{\bf{G}_2}\ra$ and the hexatic order parameter $\la \psi_6 \ra$ for the equilibrium phase transition and the laser induced freezing transition. The green dashed line denotes solid melting point $\la \psi_{{\bf G}_2}^m \ra=0.376$, and the black dash-dotted line denotes the hexatic melting point $\la \psi_6 \ra=0.06$.[]{data-label="fig:op_free_lif"}](ops_appen2_plot.pdf){width="50.00000%"}

In this section, we show the results for the relaxation time-scales of the solid and the hexatic order parameters after withdrawing an external potential commensurate with the system density under which the system is initially equilibrated. This time-scale at different densities are determined from separate simulations. The initial equilibration is performed under a time-independent external potential of the form given in \[sec:model\_simulation\] with $V(t)=\epsilon$ over $10^7$ simulation steps. Thereafter, the external potential is removed and the evolution of the solid and the hexatic order parameters are measured over time.

In \[fig:relaxation\_frequency\] we show the time evolution of these quantities at mean densities $\rho \sigma^2=1.04$ (figure $a$) and $0.98$ (figure $b$). At densities higher than the equilibrium melting point $\r_m \s^2 \approx 1.01$, the solid and the hexatic order parameters, $ \psi_{\G_2}$ and $\psi_6$, decay to finite values, indicating order even in the absence of external potential (\[fig:relaxation\_frequency\]($a$)). In contrast, at densities $\r < \r_m$, the solid and the hexatic order parameters vanish with time (\[fig:relaxation\_frequency\]($b$)). The time-scale of such decay, indicated in \[fig:relaxation\_frequency\], gives the estimate of the relevant relaxation time. If a time-independent external potential switches with a rate slower than this time-scale, the system will have enough time to [*equilibrate*]{} to instantaneous potential profiles. The relaxation frequency $f_s$ is inverse of this time-scale, and has been indicated by the dash-dotted line through $\blacksquare$ symbols in the phase diagram \[fig:phase\_diagram\].

Equilibrium Phase Transition {#appendix:PT}
============================

One way we identified the solid melting in the main text is by choosing an appropriate cutoff value of the solid order parameter $\la \psi_{\G_2} \ra$. For this, we used the equilibrium melting point in the absence of external potential. To demonstrate this, we perform separate molecular dynamics simulations of particles interacting via the soft-core potential. The temperature of the system is kept fixed at $T=1.0 \e/\kb$ using a Langevin heat bath. The corresponding variation of $\la \psi_{\G_2}\ra$ and $\la \psi_6 \ra$ with the mean density is shown in \[fig:op\_free\_lif\]($a$). The order parameters change continuously from low to high values. The solid melting point is separately determined by following the pressure-density curve and the density correlation function as in Ref. (data not shown). This melting point is found at density $\rho_m \sigma^2 \approx 1.014$, where the solid order parameter $\la \psi_{\G_2}^m\ra \approx 0.376$. The hexatic order remains significantly large even at lower densities. The melting of hexatic to liquid is identified from the change in correlation of the hexatic order parameter $\la g_6(r) \ra$, which transforms from an algebraic decay in the hexatic phase to an exponential decay in the liquid. The hexatic- melting is obtained at $\rho \sigma^2 \approx 0.998$, where $\la \psi_6 \ra \approx 0.06$. In \[fig:op\_free\_lif\]($b$) we show a similar plot for the two order parameters but in the presence of a time- independent potential of the form $U_{ext}(x,y)$ given in \[sec:model\_simulation\] with $V(t)=\epsilon$. In the presence of this potential, the hexatic does not melt in the regime of $\r \s^2 \geq 0.94$. The external potential maintains a significant hexatic order, with a value greater than $\la \psi_6 \ra =0.51$, although the solid order does drop below $0.376$ at $\r \s^2 \approx 0.96$.

![ Pair correlation functions in the driven 2D colloidal suspension at densities $\rho\sigma^2=1.04$ in the top panels ($a$-$c$) and $\rho\sigma^2$=0.98 in the bottom panels ($d$-$f$) corresponding to three different frequencies $f\tau=0.01$ (left column), $f\tau=1$ (middle column) and $f\tau=30$ (right column). []{data-label="fig:gxy_plot"}](gxy_combi_1.png){width="0.95\linewidth"}

Non-equilibrium melting and pair correlation {#appendix:pair_correlation_function}
============================================

The two dimensional pair correlation functions $g(x,y) = \la \r(x,y) \r(0,0)\ra/ \la \r \ra^2$ at $\rho \sigma^2=0.98$ and $\rho \sigma^2=1.04$ are shown in \[fig:gxy\_plot\] at the three representative frequencies: low ($f\tau=0.01$), intermediate ($f\tau=1$) and high frequency ($f\tau=30$). The figures show the correlations over a length scale of $\pm 9\s$. While the clear contrast in \[fig:gxy\_plot\]($a$), ($c$) and ($f$) demonstrate the triangular lattice symmetry, the local diffused approximately triangular structures in the other $g(x,y)$ figures are characteristic of the hexatic phase.

![image](gx_near_PB5.pdf){width="\textwidth"}

The component of the pair correlation function along the minima of the potential $g(x,0)-g(\infty,0)$ changes from a power-law to exponential decay with changing ratcheting frequency, identifying the melting point of a quasi- long ranged ordered solid. In \[fig:gx\_trans\], $(a)$–$(h)$ show the high-frequency melting, while $(i)$–$(p)$ show possible melting at low frequencies.

At $\r \s^2 \lesssim 1.03$ and low frequencies, the decay of the pair correlation function remains always exponential, identifying an absence of transition. A crossover to an algebraic decay in $g(x,0)-g(x,\infty)$ appears at $\r \s^2 \gtrsim 1.03$, resulting in solid- hexatic transition points. The phase boundaries displayed in \[fig:phase\_diagram\] ($a$) are consistent with the transition points obtained from this analysis.

@ifundefined

[39]{} \[1\][\#1]{} \[1\] \[2\] \[3\] \[3\]

F. Julicher, A. Ajdari and J. Prost, *Reviews of Modern Physics*, 1997, **69**, 1269–1282 R. D. Astumian and P. Hänggi, *Physics Today*, 2002, **55**, 33 P. Hänggi, *Reviews of Modern Physics*, 2009, **81**, 387–442 P. Reimann, *Physics Reports*, 2002, **361**, 57–265 P. Brouwer, *Phys. Rev. B*, 1998, **58**, R10135–R10138 R. Citro, N. Andrei and Q. Niu, *Phys. Rev. B*, 2003, **68**, 165312 K. Jain, R. Marathe, A. Chaudhuri and A. Dhar, *Phys. Rev. Lett.*, 2007, **99**, 190601 D. Chaudhuri and A. Dhar, *EPL (Europhysics Letters)*, 2011, **94**, 30006 D. Chaudhuri, A. Raju and A. Dhar, *Phys. Rev. E*, 2015, **91**, 050103 D. Chaudhuri, *J. Phys. Conf. Ser.*, 2015, **638**, 012011 D. C. Gadsby, A. Takeuchi, P. Artigas and N. Reyes, *Philos. Trans. R. Soc. B Biol. Sci.*, 2009, **364**, 229–238 L. Faucheux, L. Bourdieu, P. Kaplan and A. Libchaber, *Physical Review Letters*, 1995, **74**, 1504–1507 B. Lopez, N. Kuwada, E. Craig, B. Long and H. Linke, *Physical Review Letters*, 2008, **101**, 220601 P. Tierno, P. Reimann, T. H. Johansen and F. Sagués, *Physical Review Letters*, 2010, **105**, 230602 P. Tierno, *Physical Review Letters*, 2012, **109**, 198304 J. Rousselet, L. Salome, A. Ajdari and J. Prost, *Nature*, 1994, **370**, 446 S. Leibler, *Nature*, 1994, **370**, 412 C. Marquet, A. Buguin, L. Talini and P. Silberzan, *Physical Review Letters*, 2002, **88**, 168301 I. Derényi and T. Vicsek, *Physical review letters*, 1995, **75**, 374 I. Derényi and A. Ajdari, *Physical Review E*, 1996, **54**, R5–R8 R. Marathe, K. Jain and A. Dhar, *J. Stat. Mech. Theory Exp.*, 2008, **2008**, P11014 S. Savel’ev, F. Marchesoni and F. Nori, *Phys. Rev. E*, 2004, **70**, 061107 A. Pototsky, A. J. Archer, M. Bestehorn, D. Merkt, S. Savel’ev and F. Marchesoni, *Phys. Rev. E*, 2010, **82**, 030401 S. Savel’ev, F. Marchesoni and F. Nori, *Phys. Rev. Lett.*, 2003, **91**, 010601 D. Chakraborty and D. Chaudhuri, *Physical Review E - Statistical, Nonlinear, and Soft Matter Physics*, 2015, **91**, 050301(R) A. Chowdhury, B. J. Ackerson and N. A. Clark, *Phys. Rev. Lett.*, 1985, **55**, 833 Q.-H. Wei, C. Bechinger, D. Rudhardt and P. Leiderer, *Physical Review Letters*, 1998, **81**, 2606–2609 E. Frey, D. R. Nelson and L. Radzihovsky, *Phys. Rev. Lett.*, 1999, **83**, 2977 D. Chaudhuri and S. Sengupta, *Physical Review E*, 2006, **73**, 11507 J. M. Kosterlitz and D. J. Thouless, *J. Phys. C*, 1973, **6**, 1181 B. I. Halperin and D. R. Nelson, *Phys. Rev. Lett.*, 1978, **41**, 121–124 A. P. Young, *Phys. Rev. B*, 1979, **19**, 1855 S. C. Kapfer and W. Krauth, *Physical Review Letters*, 2015, **114**, 035702–5 D. Frenkel and B. Smit, *[Understanding molecular simulation: from algorithms to applications]{}*, Academic press, NY, 2002 J. Lahtinen, T. Hjelt, T. Ala-Nissila and Z. Chvoj, *Physical Review E*, 2001, **64**, 021204 E. Falck and J. Lahtinen, *The European Physical Journal E*, 2004, **13**, 267 P. M. Chaikin and T. C. Lubensky, *Principles of Condensed Matter Physics*, Cambridge University Press, Cambridge, 2012 W. Mickel, S. C. Kapfer, G. E. Schr[ö]{}der-Turk and K. Mecke, *J. Chem. Phys.*, 2013, **138**, 044501 W. Qi, A. P. Gantapara and M. Dijkstra, *Soft Matter*, 2014, **10**, 5449–5457
