Self-organising phenomena in 2D complex plasma simulations withnon-mono dispersed dust size distributions

Complex plasma with a variety of continuous and discrete dust grain size distributions are simulated in 2D with molecular dynamics simulations with radial geometry to determine differences in self-organizing phenomena to more realistically represent the actual in situ variations in dust-size. The standard deviation of particle size σ(a) strongly correlates with phase separation and coupling parameter Γ for all distribution types. We observe local differences in bond order parameters and Voronoi diagrams for different size distributions, and our results suggest that phase transition is affected by continuous size distributions, particularly in the binary distribution case. Simulations with discrete size result in artifacts and discontinuities that are not found in the continuous distributions. The use of continuous distributions is observed to be beneficial both for more realistic approximation of complex plasma experiments and to study systems of strongly coupled particles in general.


Introduction
Self organizing phenomena govern the spontaneous global order in chaotic systems which determine the nature of evolving structure, or lack thereof.Phase transitions and phase separation have been studied in coupled systems of charged particles in the fields of condensed matter, statistical mechanics, plasma and complex plasma physics [1][2][3][4][5][6][7][8].Strongly coupled systems of charged particles, atoms, molecules and dust operate at different spatial and time scales, however, the self-organizing mechanisms are the same.Dusty plasma experiments allows us to study such strongly coupled systems at the scale of millimeters and seconds [9], which has significant practical advantages.This scale-invariance also applies for simulations, which is obvious once units are normalized, with the benefit of superior diagnostics.Hence, dusty plasma experiments and simulations are a great tool for studying self-organizing in strongly coupled systems of charged particles across many disciplines [1,4,7,[10][11][12].
A binary dusty plasma consists of neutral gas, electrons, ions and two species of dust with different radii a 1 and a 2 respectively, which is characterized by the size disparity a a a ) ¯where a 2 > a 1 .In said studies of binary dusty plasma, size disparity has shown a linear relationship with phase separation in both simulations and experiments [23,24].However, the degree of phase separation was different for the same particle radii in in situ experiment and simulation.It is non-trivial to achieve perfect agreement with experiment and simulation as simulations inherently rely on assumptions and simplifications.For instance, dust and plasma dynamics exist in different time-scales and can not be easily simulated at the same time, and as such the simulations in [23] and our simulations assume a static plasma unaffected by the dust charging.Further, it is convention to assume that the size distribution is entirely discrete in simulations [23,25,26] whereas the dust particles used in the experiment has an error of 2% resulting in a continuous size distribution with 2 Gaussian peaks at the species' respective mean radii.In this paper we challenge this convention by investigating the effects of continuous size distributions.
It is unclear how a discrete versus continuous size distribution affects organizing phenomena in systems of charged particles.As all real dusty plasma experiments have a continuous size distribution due to particle fabrication error margins, it is crucial to understand how this error margin affects self-organizing phenomena.This is especially important if the continuous size distributed dusty plasma experiment is used to infer a general understanding of strongly coupled systems of charged particles with naturally discrete charge, mass and size (such as an electron).
To study the effect of a continuous size distributions we simulated discrete binary and various continuous size distributions for dusty plasma with a simplified radial 2D geometry with a newly developed code PPDyn.The set of size distributions can be found in table 1 and visualized in figure 1.
Further, we developed a new way to numerically quantify phase separation for continuous size distributions, extending the methods for binary species presented in [23].To compare our results with systems of strongly coupled charged particles in general, we have used typical diagnostics such as coupling parameter, orientational bond order, translational bond order, voronoi diagrams as well as a per particle force contribution breakdown.

Simulation setup and plasma forces
The dusty plasma is simulated in two steps, where first the plasma was simulated by [23] without any dust in a plasma chamber modelling the geometry and conditions of the real experiment in [23,24].This simulation was done with the SIGLO code (SIGLO-2D version 1.1, Kinema Software 1996-2003) [27], which adopts a fluid approach solving the hydrodynamic transport equations for charged particles.
Next, the dust particles are simulated in a static plasma, where the plasma induced forces are collected from the first simulation.Binary dusty plasmas have been simulated with molecular dynamics codes [23,25,28,29].However, we have developed a new fully parallel and C-compiled Python based code called PPDyn tailored specifically to the case of dusty plasmas.Specifically, our codes provide easy implementation of different size distribitions such as: • Discrete where all the particles are defined exactly by a given set.This is accurate for particles of discrete size and charge.
• Gaussian distribution about one or multiple mean radii, which emulates manufacturing errors in real dust experiments such as [6,23,24,27].
• Uniform distribution within a range, which is similar to space plasmas such as Saturns rings [30] In this investigation we have simulated 6000 particles for one second with 3000 timesteps in 2D with particle size distributions according to table 1.The code uses the velocity Verlet method [31] to move the particles.The simulated chamber has a simple radial geometry extending 15mm radially from the center of a circular 2D chamber.The geometry in the real experiment in [23] has a wall at about 18 mm, however the potential well, as seen in figure 2 prohibits dust particles from exceeding past 12 mm radially.As such, the edge effects at the physical walls of the chamber are neglected, as the Debye length is 300μm.The dominating forces considered in the simulations are the ion drag force, the electric field, neutral drag, and inter-particle repulsion.The plasma parameters can be found in table 2.

Implementation of forces
The plasma is not directly simulated in PPDyn and the forces acting on the dust particle have been adapted from previous work [23] as in figure 2. The forces are spatially interpolated to the particle position.The data describing the forces impacting the dust are collected from SIGLO by [23], where plasma has been simulated in a chamber with the same neutral gas density and geometric dimensions as used for the simulations in this investigation.The plasma interaction with the dust particle of radius a is calculated according to the ion drag F ion and the electric field in the chamber F EL from [23] for a = 3.5 μm in figure 2(a) and a range of radii in (b).
As seen in figure 2(b) the net force F tot reaches an equilibrium around r = 11 mm depending on the particle radius.To save computational resources the dust particles are loaded at rest close to r = 11 mm as seen in figure 3a.

Dust charge
The dust charge is non-trivial to solve analytically, and is perhaps best solved with a combination of experiments and numerical methods.Regardless, we will provide a brief explanation of dust charge, and the most important factors impacting it.The ion and electron currents I i and I e with respective number densities n e and n i to a particle with radius a and charge f p < 0 in a stationary Maxwellian plasma is given by the well known Orbital Motion Limited (OML) model For Argon, and T e /T i = 100 we obtain 2.41 f = -ˆ, such that our OML charge is Z d ≈ 6300.However, we do not have a stationary plasma, and as such ion streaming should be taken into account.For v v 0.5 200 it has been shown that Z d increases as much as by a factor of 2 in [33].We have v v , meaning that a greater Z d due to streaming ions is appropriate in our situation.From [34] ion-neutral collisions will modify the ion curernt Where ℓ mfp is the mean free path for ion-neutral collisions.As our particle is negatively charged, the resulting ion current due to ion-neutral collision will be larger than in the collisionless case, overall reducing our charge number Z d .Lastly, we consider the dusty plasma effect, or  In (a) we plotted the forces felt by a dust particle of size a = 3.5μm from ion drag, F ion , and the electric field, F el , as a function of distance from the center of the chamber r.The values are from the SIGLO simulation in [23].Reprinted (figure) with permission from [23], Copyright (2021) by the American Physical Society.Here, the force from the ion drag is uniformly increased by a factor of 2 in the implementation (but not in the plot), as also seen in the total force being a sum of the forces for F tot .Note that r is considered in the radial direction, effectively making this a 1D case.In (b) we have plotted F el + F ion for a range of particle radii to show their different y-intercept (where F el + F ion = 0).
Maxwellian background plasma, and semi-infinite cloud-size [35], and cloud thickness H where H > > λ S which is not our case.In fact, we have plasma production in the middle of our cloud, and it is unlikely that our 3mm thick cloud can entirely contain the plasma discharge such that no electrons enter the cloud due to the dusty plasma effect.We suggest that electron screening from the dust cloud will push against the electron discharge in the center of the chamber, such that the cloud is pushed further out radially, however it will not substantially alter the results.Further, [23,24] have good agreement between simulation and experiment without accounting for this effect in their MD simulation.
Overall, it is hard to know the exact realistic charge, as many effects and parameters affect not only the dust charge, but also each other.With the effects mentioned, and considering that we are basing our charge of another simulation which already has its own set of assumption, we suggest that the actual charge is probably within an order of magnitude of the OML value of Z d = 6300.As such, we started out with the results used in [23], and did some changes to optimize for computational time.In [23] for a = 3.5μm the charge number was found to be Z d = 9460 for the SIGLO simulation.However, using [23] as our benchmark, and to lower computational expense, we increased the dust charge by a factor of 5. Further, increasing Z d did not impact the phenomena of interest in this study, as long as the relative charges remain the same, however it substantially cut down on the computational time needed to reach a steady state.Additionally, slightly artificially increasing the neutral drag caused the system to reach the steady state faster without impacting results.
One major limitation of our dust charge model is that the dust charge is static.Realistically, dust is charged dynamically according to the plasma conditions at the particle position.From equations (1), ( 4) and (5) we see a dependence on temperature, number density and mass for each species, in addition to the individual velocity distributions.Future work should consider a dynamic charge model, possibly implementing a hybrid method with a simple fluid model of plasma and dust as point particles.The effects of charge fluctuations are found to cause dissipative and instability mechanisms for ion waves in the plasma [37,38].However, we assume that the greatest contribution to our dust dynamics is the size differences' effect on Z d .Even if we consider a dynamically charging model with dependence on r in our radial symmetry, a binary species will separate according to their size, as the larger particle always will obtain more charge.
is the electric field in the simulated plasma chamber.The force acting on the dust particle due to the electric field is given as [23] Table 2. Plasma parameters for the SIGLO simulation in [23] and our dusty plasma simulation.The electron densities are approximations for the simulation box, where the detailed plot can be found in [23].r min is the distance to the closest neighbor.T e and n e are electron temperature and number density, Ti is ion temperature, P is neutral gas pressure, λS is the screening length, Urf is the electric potential difference and ν rf is the rf frequency and ℓ mfp is the mean free path of ion-neutral collisions.a ¯is the mean dust particle radii, n d is the dust number density, Z d is the dust charge of particles with mean radius a ¯, and r min á ñis the average distance to the closest neighbouring dust particle.
Plasma parameters The collection force is [23] F amvnv e m v 1 2 10 Where a is the dust particle radius, m, v and n is the mass, velocity and densities of the respective species where the subscript i is for the ions and In our case, Hutchinsonʼs approximation T e = 100T i fits well [23] for electron and ion temperatures respectively, and the collisionless orbit force is [39] F anGu e m v 8 l n 1 3 where u = v i /v T,i [23].The Coulomb logarithm lnL is given by where λ eff is the effective screening lentgth taken from [39] such that k T m v a 1 2 17 .Finally, the total force found in the plasma chamber simulations impacting the dust particles of radius a = 3.5μm as seen in figure 2 is given by F total = F el + 2F ion , where F ion is from SIGLO simulation in [23] seen in figure 2.

particle-particle repulsion
The particle-particle repulsion between two dust particles at distance r is governed by the well known screened Coulomb potential Where Q is the charge of the respective particles and λ S = 300μm in all the simulations [23].From table 2 the screened Coloumb potential is appropriate.The particle repulsion is calculated between all particles.

Neutral drag
The drag force due to neutral gas is is based on the Epstein formula [44] N v m 1. 33 19 Where N n is the number of neutrals given by PV = NRT, m n is the mass of the neutrals and The force felt by the dust particle is given as The forces are first computed for a particle with radius a = 3.5μ m, then scaled according to each force's dependence on a.

Results
Here we present a variety of diagnostics to quantify the different size distributions' phase separation, phase transition, structure, coupling parameter and order parameter.Comparisons to other works are presented alongside the results, while a broader discussion is found in the discussion section.

Phase separation
First, we investigated the particle distributions' effect on phase separation.Phase separation has been shown in binary dusty plasma experiments in micro-gravity as well as simulations [23,24], where simulations using discrete binary size distributions had less separation than the experiments.Here, continuous and discrete size distributions have been simulated to investigate this discrepancy further.¯¯is plotted for different size disparities ò where a a a 22 In our simulations a 2 and a 1 are the mean radii of species 1 and 2 and a ¯is the mean value of a 2 and a 1 .We clearly see a faster and greater degree of phase separation for larger size disparity in figure 4, which corresponds to results in previous work [23,24].

Continuous distribution
To quantify the order of separation in the continuous case, a method different than b b ¯¯is used as essentially, 6000 species with 6000 different radii were being dealt with.In figure 5 the particle size is plotted versus radial position for the continuous distributions' smallest (a), (c), (d) and largest (b), (d), (f) size disparities.In figure 5 we observe that for the largest size disparities particle size seems approximately proportionally correlated to radial position indicating phase separation.However, for the smaller size disparities, it is not clear whether this is the case.
We also observe a larger dispersion close to r = 12 mm in figures 5(b) and (f).As the dispersion is larger at the extremes of the Gaussian curve, we should expect more dispersion at the edges for a conmpletely separated Gaussian particle size distribution.Additionally, the outermost shell needs the most particles to fill, and again there being fewer and more dispersed particles at the extremes of the Gaussian distribution contributes to this effect.Lastly, a sharper field gradient f   at the edges of the potential well should also contribute somewhat to this effect, and largely explains why we don't observe the same at the inner edge as there the gradient is very low.However, we don't see more dispersion at the outermost edge in the uniform distribition in figure 5(d), so the gradient likely only contributes to less dispersion at the inner edge.
In order to quantify the phase separation also for continuous size distributions we divide space radially into 20 segments r i .Then, the average particle radii a i ¯is calculated for the corresponding segment radially.Now we can find Where C i is some constant.With the discrete particle distributions we have used ò as in equation ( 22) as a measure of disparity.However, for continuous distributions we instead compute the standard deviation of particle size σ(a) for each size distribution as a measure of size disparity.
In figure 6 we have plotted β for the different size distributions, and it suggests that the standard deviation σ (a) predicts β regardless of distribution type.Macroscopically this means that the degree of phase separation is independent of distribution type, and only dependent on overall standard deviation of the particle sizes.The deviations from a perfect linear relationship has been likely due to how the forces scale non-linearly with a, and each distribution type caused a different local distribution during and towards the end of the phase separation.
Further, the separation forces are a result of F ion , F el and interparticle repulsion, which for each particle depends on it's own size, and the size of neighbors.As such, we suggest heterogeneity that depend on the size distrbutions' individual topoligy's effect on how particles populate the space with respect to the static forces.In figure 5 we see how the particle size distribitions' topoligy manifests itself in the radial position.Hence, the linear fit of β is weighted according to N d , which is the number of particles at a small section radially with width δr, to take some of the size distribitions' different radial position into account.Further, we suggest that some particle size distribution topology can be incorporated in an improved model for β, but as its effect is not dominating developing such a model is outside of the scope of this paper.However from figure 6 the variation in β is not dominated by this heterogeneity, and as such the weights in our fit is considered sufficient to show σ(a) as a predictor of separation regardless of size distribution.¯¯over the total time B is shown for different ò.We observe B ∝ ò which agrees with previous work [23,24].
For instance, in distributions with a higher proportion of larger particles, the outermost particles have a greater repulsion, effectively pushing the rest of the particles further towards the center and offsetting the position of the static forces seen in figure 2. This effect is larger for the uniform and binary Gaussian distribution, and smaller for the Gaussian.
Micro-gravity experiments have been conducted with binary dusty plasma with size disparities ò = 0.05 [45] with great degree of separation.In figure 3 we have ò = 0.09 and comparable degree of separation, at least from what is discernible visually in the two cases.
Notice that the binary Gaussian distributions have a greater σ(a) than its binary equivalent in figure 6.Here, equivalent refers the discrete size disparity, where the separation between the binary species' radii are equivalent to the separation between the mean values for the binary Gaussian distributions.In other words, the real experiment with some error for the binary species' radii, will for the same discrete size disparity have a greater σ (a), resulting in more separation.As such, one should expect more separation in a real experiment than in a simulation with discrete sizes for the same discrete size disparity.
In [23] the experiments also observed larger degrees of separation than the simulations, possible due to the same effect.This effect is due to half of each species having a larger size disparity due to the two Gaussian Figure 5.The radial distance from the center of the chamber r versus particle radius a for the smallest 'min' and largest 'max' size disparity for binary Gaussian, uniform and Gaussian distributions.For the smaller size disparities we observe bonds forming at the outermost r.For the largest size disparity we see a correlation between a and r, where the uniform distribution suggests a ∝ r. ¯is the average particle radii within one of 20 segment of radial distance from the center of the chamber r i where r is adapted to contain every single particle such that r r r min , max i j j Î ( ( ) ( )) versus the standard deviation of the particle radius σ(a).The line was fitted with weights proportional to N i where N i is the number of particles in theith segment.We observe β ∝ σ(a), which agrees well with B ∝ ò for the discrete binary case.distributions, and the other half having a smaller size disparity.From our results it appears that this does in fact not cancel out, and overall leads to larger degree of separation.It is unclear whether the discrepancies between simulations and experiments can be entirely explained by discrete vs continuous size distributions.All dustplasma effects are neglected due to different scales in the simulations, which likely is a considerable contribution to the discrepancy.

Coupling parameter
Systems of charged particles are considered coupled when the Coulomb energy per particle is greater than the thermal energy.The respective energies' ratio is called the coupling parameter Γ and for a screened Coulumb potential is given [46] Z e k Tr exp 25 Where T is the temperature of the dust particles, 〈Z d 〉 is the average charge number and r WS is the 2D Wigner-Seitz radius given by Where η is the number density, which for a 2D system is given per area, and κ = r WS /λ S where λ S is the screening length.A system is strongly coupled for Γ > 1, and crystalline structures are found in systems of charged particles with Γ > 168 ± 2 [2].In table 3 the coupling parameter is computed for the different distributions in our simulations, where larger σ(a) correlates with lower coupling parameter.
From figure 7, which shows the coupling parameter Γ as a function of σ(a), we observe that a log 1 s G µ ( ) ( ).As such, increasing the standard deviation in particle size exponentially degrades coupling.For the points in figure 7 at σ(a) ≈ 0.2μm there is large variety in the coupling parameter, suggesting that the different distributions have some effect on how Γ depends on σ(a).Further, linear proportionality agrees better within each distribution type separately than as an overall trend across all distributions.Overall the uniform distribution has the largest Γ and the smallest σ(a).
From equation (25) the dominating variable is the temperature of the dust particles, which interestingly infers that T a log s µ ( ( )).If we consider temperature as a measure of entropy, and large σ(a) as a disruption in the order, it can help understand this relation.One could argue that the differences in temperature is a result of the simulation not reaching a steady state.However, the runs are initialized the same way and are run for the same amount of time, with the only difference being the size distribution.As such, the relation between σ(a) and temperature must be the result of non-linear effects from the sum of particle interactions.
However, we here assume constant particle charge, where realistically the charge is dynamic with fluctuations as discussed in section 2.1.1.We suggest that charge fluctuations will induce random motion in the dust species, effectively increasing the dust temperature reducing Γ.Further, the dynamic charging is strongly space dependent, where local variations in plasma conditions will determine the dust charge.Charge fluctuations might perturb the structure locally or globally, inducing dust acoustic waves, where dynamic charge has proven dissipative in [37,38].
For σ(a) < 0.1 in figure 7, which corresponds to a ratio of about a a 0.03 s » ( ) ¯, we observe Γ corresponding to a solid phase, indicating a solid Coulomb crystal [8].However, for σ(a) > 0.1 the coupling parameter indicates a liquid phase.If we compare the discrete binary and the binary Gaussian distributions, the relation between σ(a) and Γ is very different, particularly close to the phase transition.For the solid phase the binary Gaussian has slightly greater Γ than the discrete binary, while for the liquid phase the binary Gaussian has significantly lower coupling parameter (by a factor of 4).The binary Gaussian distribution has the largest difference in coupling parameter from solid to liquid (column (a) and (b) in table 3).We argue that the binary Gaussian distribution is the most complex, and as such can be expected to be the least predictable.As such, phase transitions should be expected to behave quite differently in simulations with discrete binary distributions compared to real experiments with two Gaussian distributions.Cylindrical dusty plasma experiments with layered ring crystals have shown coupling parameters in the range 370-490 for κ in the range 1.7 to 2.4 [25].Our simulations have κ = 0.4, where Γ/κ diagrams grow exponentially in [8], resulting in good agreement between our Gaussian distribution with error of 2% and the results in [25].However, the size error of the dust particles were not shown in [25], only the mean radii a m 2.55m = ¯.Numerical studies have shown a exponential proportionality of both κ and Γ to a/λ S [26], suggesting that continuous size distributions have non-linear effects on particle coupling.Our results suggest a similar effect as seen in the differences already discussed in figure 7 and table 3, however the studies in [26] did not use a continuous size distribution but particles with radius a = 1μm, which is discrete and smaller than our particles.Overall, we show that the distribution type impacts coupling and phase transitions in a non-linear way, which has been inferred from previous works [8,26] but not explicitly shown.

Local structures
Next, the local structure is quantified and examined.In figure 5 there are bonds forming for the smaller size disparities in (a), (c), and (e), which are not present for the largest size disparities in (b), (d) and (f).This suggest that there are different local structures both radially and for the different size disparities.
Pair correlation functions, orientational bond-order plots, and Voronoi diagrams are produced to investigate the distribution types effect on local structures.The pair correlation and bond order has been used to determine the phase of strongly coupled systems [4,47].
For the translational pair correlation and orientational order plots the particles have been divided into 4 regions radially in order to compare and isolate the local differences.Due to different levels of separation, and different size distributions, different particle sizes populate each of said 4 regions.There is some overlap between the regions to compensate for low particle number in the innermost and outermost region.
It's been shown that 2D systems always have long-range fluctuations disrupting order [4,13], such that our results will not be the same as in a 3D systems.It has also been suggested that freezing is a two-stage process in 2D, first undergoing a hexatic phase, which is not present in 3D [4], which will be discussed later.

Translational pair correlation function
For a given particle it has been investigated how many particles δN i are found at a distance (r i , r i + δr) such that g r N r r Where ρ is the particle density.Averaging over all particles we get the translational correlation function g(r), which indicates the range of the order of the system and gives the average particle distance.
In figure 8 the translational correlation suggests strong translational order for smaller size disparities with at least 3 distinct peaks.For the larger size disparities the translational bond order is weaker, with at most 2 distinct peaks.Additionally, the mean particle distance is smaller towards the outside of the circle which is where the forces acting on the particles are strongest.

Orientational order
To further quantify the order we have used the orientaional bond order [48] by assuming a hexagonal structure.A perfect crystal lattice will be equiangular hexagonal in 2D [49], whereas in 3D it will have BCC (Body Centered Cubic) crystal structure [50].Hence, by quantifying how hexagonal our structure is, we also quantify the extent of a 2D crystal lattice structure.Next, we draw the bonds from a given particle to its closest 6 neighbors, where the angle θ i is the angle between each bond and some common vector pointing radially from the given particle to the center of the chamber.The orientational bond order is given by Where M = 6 and Ψ is a complex number of magnitude |Ψ| 1.Further, we evaluate the orientational alignment of the closest six neighbors for any particle pair i and k separated with distance r r r

Averaging over all particles we obtain [51]
G r r r 29 The orientational order in figure 9 has the same overall trend as the translational order; more ordered for smaller size disparity.

Voronoi diagrams
A third way to evaluate the order in the hexagon is done by quantifying how much each bond to the closest 6 neighbours with length l i deviates from the average bond length μ l .A perfect hexagon will have all l i = μ l .Hence, a bond length parameter L is given by Where M = 6 corresponding to the closest 6 neighbours.Combined with Ψ 6 , we obtain both a measure of equal bond length and equal bond angle.In figure 10 the product L • Ψ has been plotted as an overall indication of the hexagonal structure per Voronoi cell.Lets first consider the discrete binary and binary Gaussian distributions, where for the smaller size disparities in figures 10(a) and (c) there are islands of hexagonal lattices separated by comparable islands of dislocations.However, for the larger size disparities we mostly observe a hexagonal lattice with few dislocations for the discrete binary in figure 10(b), particularly along the separation border between the two species at about r = 11 mm.However, the binary Gaussian in (d) has some small islands of hexagonal shape with mostly non-lattice structuring.
Next, for the uniform distribution in figure 10(e) we observe the most structured hexagonal lattice, with only some dislocations.Unlike the discrete binary, there is no separation border between the two discrete species causing dislocations for the uniform distribution.In (f), we observe the smallest degree of hexagonal lattice structure.In both 10(e) and (f), the structure and lack thereof, respectively, appears uniform in the interval r = (10, 12) mm.
Lastly, the Gaussian distribution in figures 10(g) and (h) has the least overall structure in their columns.Further, it appears that the more structured regions in the Gaussian distributions are populating the middle of the particles' radial domain.
Overall, we observe that the different particle size distributions' topology correlate with where we find the hexagonal lattice structure radially.The structured regions are where the neighbouring particles are closest in size, which will depend on the particle displacement due to separation and the probability that the neighbouring particles will be close in size.Now, to fill the outermost layer we need the largest number of particles, which will also be the largest particles due to the separation.However, given a Gaussian particle size distribution, the size variation will be greatest at the extremes, such that the outermost layer should have less structure, as seen in figures 10(g) and (h).We observe equivalent effects for the other distributions.
Further, we observed different behavior for the discrete binary and continuous distributions.In particular, the discrete binary distribution appears to get more structured for larger size disparity while the continuous distributions suggest less structure for larger size disparities.For large size disparities the discrete binary dust is almost entirely separated, effectively resulting in two mono-dispersed species separated with a mixing layer with some dislocations as seen in figure 10(b).However, we don't see the same behaviour for any of the continuous distributions despite comparable levels of separation.
In figure 11 the average of the product 〈L • Ψ 6 〉 is plotted against the standard deviation of particle radii σ(a) for each distribution type.Here, the distribution type does not seem to matter much, leaving the parameter σ(a) as the dominating predictor of 〈L • Ψ 6 〉.Comparing figure 11 to figure 7, the coupling parameter has much more variation depending on distribution type, indicating that the coupling parameter is more dependent on each distribition's individual topological differences.
The number of higher peaks are observed to be declining with respect to r in figure 9 for all distribitions.However, the uniform and Gaussian distributions appear to retain more peaks better, particularly in the second column in figure 9.
Comparing the discrete binary and binary Gaussian distributions, we see large differences in figure 10.Futher, while the differences are smaller, we still note some differences in figures 9 and 8, such as the shape of the outermost segment (blue line) in column 3 and 4. As such, we have demonstrated that a real experiment with two Gaussian size distributions and a simulation of two discrete species can behave quite differently.This suggests that using binary Gaussian distributions in simulations can help better approximate real experiment behaviors such as in [23].On the left there are several peaks indicating a solid phase and long range order.Note the slight differences in binary and binary Gaussian, particularly for the blue line, indicating that discrete and continuous binary distributions do not behave exactly the same.

Force per particle
As outlined in section 2.1 the particles will feel the effect of the ion drag, electric field, inter-particle repulsion and neutral drag.For a system at steady state, the sum of the forces will be zero, with zero neutral drag, such that Where F repc is the repulsive force radially.Further, when F el = F ion , which is the y-intercept in figure 2, we also need the repulsive force to be zero for net zero forces at the steady state.Moving away radially from the point where F el = F ion , it follows that F repc ∝ F ion + F el assuming a steady state.
For mono-dispersed particles, the point where F el = F ion will be the same for every particle as the forces depend on particle radii a.However, for binary particles, there will be two such points, and for continuous size distributions every particle will have a unique profile for F el and F ion (as seen in figure 2).
In figure 12 the forces per particle have been plotted for the largest size disparity for the discrete and continuous binary distribution in (a)-(b) and (c)-(d) respectively.First, we notice a discontinuity for the discrete binary in (a) and (b) at the separation border slightly above r ≈ 10.5.However, for the continuous binary Gaussian in (c) and (d) no such discontinuity was found.In figure 10(d) we noted a breakup of hexagonal structure at said discontinuity in figures 12(a) and (b) which was not present in the continuous distributions.Essentially, the binary Gaussian size distribution results in small differences between the forces acting on the particles, which in any real experiment would be the case due to small fluctuations in F ion and F el as well as dynamic charging of the dust particles.As such, the discontinuity and breakup of structure at the discontinuity has entirely been a simulation artifact and not realistic.
Next, we observe that the repulsive force F repc is proportional to the sum of the static forces F ion and F el in figures 12(b) and (d).However, they're not equal which indicates that the system is not yet at rest.Further, the point where F el = F ion (sum of forces is zero), which corresponds to the discontinuities in figure 12, is clearly On the left there are several peaks indicating a solid phase and long range order.Note the differences in binary and binary Gaussian, indicating that discrete and continuous binary distributions do not behave exactly the same.defined for the discrete binary in figures 12(a) and (b) as they have the exact same size and thus the same point where F el = F ion .However, for the continuous distribution in (c) and (d) where every particle has a different point for F el = F ion the zero point has dispersion correlated with particle size dispersion.As particles will be pushed towards the F el = F ion point, and this points being different for every particle, it will result in a shear force disrupting the lattice or hexagonal structure.Ultimately, the size in our case is regulating the dust particle's charge, so F el = F ion depends on particle charge only.In a real experiment, the particles will charge dynamically, and as such the F el = F ion point will change dynamically for every particle (not to mention that F el and F ion will fluctuate), perhaps wiggling about a time averaged point.As such, the case in figures 12(a) and (b) where F el = F ion point is clearly defined is not realistic and also ultimately a simulation artifact.

Discussion
In our diagnostics we observe both similar and different behaviour for discrete and continuous size distributions.We have already discussed and made some comparisons with other work for each specific diagnostics type.Here, we will take a holistic look at our results and discuss the impact of continuous vs discrete size distributions in the context of existing simulations and experiments in strongly coupled systems.
Figure 10.Combined L • Ψ 6 calculated for each particle to mapped to color each Voronoi cell for the smallest (a), (c), (e), (g) and largest (b), (d), (f), (h) size disparity for each of the size distribution types.The cropped region is zoomed in to the left part of the particle ring.For the continuous size distributions we have higher L • Ψ 6 for the smallest size disparity.However, for the discrete binary distribution this trend is reversed with higher L • Ψ 6 for larger size disparity.

Hexatic phase
While the system is not a crystalline solid, not an isotropic liquid, there exists in 2D a hexatic phase in between according to KTHNY theory [4,13,14,16].The requirements for the hexatic phase are conventionally short range translational order where g(r) ∝ e − r/ ξ and quasi-long orientational order where G 6 (r) ∝ r − ν with 0 < ν 1/4 with some dislocations disturbing the hexagonal structure [47,52].Figures 8 and 9 suggest the hexatic phase for the smallest σ(a) at the outermost region r = 11.5mmparticularly for the uniform distribution seen in figure 8(e) and 9(e).We were able to fit approximate lines with ξ = 0.2 and ν = 2. Here, ν is outisde of its acceptable range, however that would require translational pair order which is impossible with our geometry.As such, we suggest that a hexatic phase from earlier is justified.
In figure 13 the distortion of hexagonal lattices are shown in a simple diagram.Now we will try to quantify these errors with some approximations.Let's assume we have a shell of evenly distributed particles, such that the number of particles in the shell N , which in our case should not be more than 1/4 of d before results of the translational pair order stops working, which is around 2 mm.The translational pair order will deterioate gradually up until this limit, such that translational pair correlation with longer range than this limit is also impossible for a perfect crystal.Further, the distortion in orientational order is given by r tan d r ), which in our case is about 0.01rad.As our orientaional order operates with multiples of π/6, this is an error of only 2%.As such, the translational pair order is much more significantly impacted by the circular distortion of the hexagonal lattice, than the orientational order.
Note that the uniform distribution has the lowest σ(a) with strongest bond order correlation, and the single Gaussian has the highest σ(a) with the weakest bond order correlation.As such, the hexatic phase might be correlated to the overall particle size distribution, where high enough size disparity σ(a) makes the hexatic phase impossible.However, when the discrete binary distribution is completely separated, we obtained two monodispersed species inhabiting different regions radially: One with only species 1, one mixing layer and one with only species 2. The two regions with only one species present will have σ(a) = 0 locally, while the overall σ (a) ≠ 0. Hence, looking only at macroscopic predictors such as σ(a) might suggest similar behaviour for discrete and continuous size distributions while in reality the local structures are different.Similarly, the local structure in figures 10(c) and (d) suggest that the discrete binary distribution is more ordered locally for higher size disparity σ(a) outside the mixing layer at r > 10.5, while the macroscopic average 〈L • Ψ 6 〉 in figure 11 suggest less order for higher σ(a).structures of the particle size distrbutions.Hence, the coupling parameter depends on the local structural differences we see in figures 5 and 10.Further, there are some issues related to the bond order parameters [53], particularly in the 6 closest neighbours.The particles at the outermost and innermost layer will not have 6 neighbours, but 4. As such, the innermost and outermost particles will falsely indicate lack of structure.Others have improved on this with a more sophisticated choice of neighbourhood [53] and weighted bond order [15].
The average particle distance to the closest neighbour depends on the repulsive force F repc ∝ F ion + F el .This means that the translational order is not expected to be the same everywhere, but will vary radially depending on F ion and F el .As the profile of F ion and F el depends on a and is unique for every particle, accounting for and correcting for this non-linear particle repulsion effect in the translational correlation is non-trivial.Further, the radial symmetry fundamentally does not work with a perfect hexagonal structure, unless the hexagon is distorted to fit a circular geometry.Regardless, this effect is also not taken into account for the bond order parameters.
Overall, the complexity of the inter-particle forces and differences in the F ion and F el profile for each particle lead to a complex particle distribution in space over long distances that escapes particularly the translational bond order (as it assumes the same inter-particle distance everywhere).As such, the bond order indication that the particle distributions have similar structure from figures 9, 8 and 11 might be subject to the common shortcomings of the order parameters outweighing their differences.

Phase separation vs bond order
From figures 11 and 6 we see that both β and 〈L • Ψ 6 〉 can be predicted by σ(a).This suggest that phase separation indicates the structure of the system.Even for the discrete binary distribution for large size disparities where the two species are almost completely separated (as indicated by the discontinuity in figure 12) we still see separation indicating less structure.
Phase separation indicates a difference in forces acting on the particles, according to size a in our case.Structure on the other hand indicates how regular the inter-particle forces are, where a perfect crystal has the exact same inter-particle forces everywhere.This suggests that any sign of partial phase separation is detrimental to a perfect crystal lattice.For particles to separate, we need that the separation force is greater than the particle coupling force keeping the particles in their place in the lattice.Hence, separation indicates that said force upholding the lattice structure is too weak.
However, it is possible to completely separate in binary species such as multi-component fluids before phase transitions.But, as shown in our results, real dust with binary Gaussian distributions do not appear completely separate and then form crystals the same way discrete binary systems do.Further, as we observe phase separation even for the single peak Gaussian with σ(a) at 2% any attempt at studying crystal-like structures should be aware that the particle size error has a significant effect on what type of structures one should expect to see.

Comparing discrete and Gaussian binary distributions
Binary dusty plasmas are typically simulated with discrete sizes [23] and compared to experimental data where sizes are continuous.If we consider the discrete binary and binary Gaussian distribution simulation runs we can evaluate some of the similarities and differences as a result of continuous vs discrete sizes.In table 4 we see that using a continuous distribution results in overall larger σ(a), something to keep in mind for metrics where σ(a) appears to be a good indicator.
In [23,24] phase separation in discrete binary simulations are compared to micro-gravity experiments where greater separation is observed in the experimental data.In our experiments, we observe a general trend of greater separation for the binary Gaussian distribution.However this is likely due to continuous and discrete distributions having different σ(a) for the same ò, as seen in table 4. Further, if σ(a) is taken into account, it appears that discrete and continuous binary distributions largely have the same degree of separation according to our metrics.
For the coupling parameter however, we see the largest diversion between the discrete and Gaussian binary distribution as seen in figure 7. Particularly, the binary Gaussian has a very large increase in Γ from the liquid to the solid phase compared to especially the discrete binary, but also all other distributions.This indicates that phase transitions may be significantly affected by continuous vs discrete implementations, something we also see in figure 10.Further, the discrete and binary Gaussian have different scaling laws for Γ and σ(a) in figure 7, suggesting that coupling, and phase transitions, work differently for Gaussian distribitions.
The reduced coupling parameter due to large σ(a) is similar to the phenomena known as 'freezing point depression' [54], where introducing impurities into a liquid (such as salt in water) effectively reduces it's freezing temperature.However, it appears that in our case the larger σ(a) actually prevents cooling rather than lowering the freezing point, as the final temperature is different for different σ(a) as indicated by Γ.Whether the mechanism is exactly the same or not is not clear, however it is clear that both continuous size distributions and impurities as described by [54] disrupts crystal formation.
Overall, phase separation seems relatively unaffected once changes in σ(a) is accounted for.However, the coupling and phase transitions appear to be largely affected by the different distributions.

Benefits of continuous vs discrete
One could argue that the point of a simulation is to investigate ideal and simplified scenarios with unrealistic assumptions for simplicity.If so, it is extremely important to be aware of such assumptions.On the contrary, the simulation artifacts due to simplification might cause unnecessary barriers when comparing simulations to real experiments.In our case, there are several benefits of the continuous vs discrete distributions 1. Particle charge changes dynamically and will typically not be homogeneous even for the exact same particle size.A particle distribution size better models dynamic particle charging.
2. Discrete binary species has force discontinuity and too well defined net zero force point as seen in figure 12 which is unrealistic.
3. A hard border between species as seen in figure 10(d) disrupts particle structure.
Overall, assuming that continuous size distributions are trivial to implement, there is no reason not to have them.It removes a lot of ambiguity especially for local structures and disregarding unrealistic simulation artifacts.
4.6.General applications of σ(a) in strongly coupled systems Our results suggest that varying σ(a) has a significant effect on phase transitions.We know that mono-dispersed Coulomb crystals undergo phase transitions according to their Γ − κ phase diagrams [8].Yet, it is unclear whether the same conditions apply for large σ(a).However, the conditions for the hexatic phase in [47,51] agree with our diagnostics for low σ(a).While we refer to σ(a) as radius, its physical effects on the system manifest themselves in the particle charge and mass' dependence on size.Further, the charge and mass affect both particle-particle and particle-field interactions.Despite strongly coupled systems of ions operating with discrete sizes, masses and charges, we suggest that fluctuations in particle-particle and particle-field interactions are analogous to σ(a).The essence of this suggestions, is that the existence of imperfections, be it in size, charge or field, will have a non-linear effect on the order of the system.Further, as our results show, statistically quantifying the imperfections allows us to predict and model their effects.

Limitations of our model
Here we will briefly discuss some limitations of our model and assumptions and possible impacts on our results.First, we assume constant charge of the dust particles.Fluctuating charge has shown dissipative [37,38], which in our case would erode structuring.Further, as we have shown that small changes in Z d due to a normal distrbution of particle size impacts self organizing phenomena and phase tranitions, we can assume that small charge fluctuations in Z d will impact the same processes.Second, we assume negligible sheath edge effects, as our dust particles are populated in a potential well situated many Debye lengths away from the chamber walls.Our experiment exists in microgravity, however other experiments with a strong potential to counteract gravity might observe significant sheath edge effects.
Third, we assume a negligible dusty plasma effect.As shown in section 2.1.1,our 2D disc of particles have a Havnes parameter P H ≈ 0.003, justifying the negligible effect.However, for a sphere of equal particle density in 3D where particle number N = 464000, the Havnes parameter is P H > > 1, meaning we have complete screening of electrons.But, our dust cloud is much smaller than what the Havnes effect assumes, and we don't have a Maxwellian homogeneous background plasma.As such, exactly how a spherical cloud is affected remains somewhat ambiguous, and should be considered future work.

Summary
We have studied the phase separation and bond order for discrete and continuous size distributions in dusty plasma simulations.Our findings can be summarized to 3 main points, highlighting differences, similarities as well as some findings independent of distrubtion type.Overall, the standard deviation of all the particle sizes seem to be a very important indicator, both for phase separation and structuring.The main findings are the following 2.There appears to be local structural differences for large size distributions and a high degree of separation for discrete and continuous size distributions.Temperature is proportional to a log s ( ( )).However, this effect is different for different size distributions at roughly the same σ(a).

5.
Phase transitions appear to be affected by discrete vs continuous size distributions, particularly in the binary case.
6. Static forces and discrete particle sizes creates artifacts in the force per particle resulting in disruptions at the species' separation border and a too stable and well defined zero force resting point not present in continuous distributions.
Further work should consider whether these findings also apply for 3D, and overcome the shortcomings of the bond order outlined in section 4.2 and [53].Additionally, perturbing the static forces F ion and F el , or dynamically charging the dust, should be compared to the effect of continuous size distributions.Continuous size distributions might satisfactory account for perturbations and dynamic charging, or perhaps other phenomena could be observed.Lastly, such effects should be compared to experimental data to determine which effects should be modelled, and which could be neglected.

Figure 1 .
Figure1.The distribution of of particle radius for the different types distributions for increasing range of particle radii a for (a), (b), (c) and (d) respectively.Blue is Gaussian, red is uniform and black is binary Gaussian.This makes up the 12 datasets used for the 2Dsimulations with distributions of radii of 6000 particles.

Figure 2 .
Figure2.In (a) we plotted the forces felt by a dust particle of size a = 3.5μm from ion drag, F ion , and the electric field, F el , as a function of distance from the center of the chamber r.The values are from the SIGLO simulation in[23].Reprinted (figure) with permission from[23], Copyright (2021) by the American Physical Society.Here, the force from the ion drag is uniformly increased by a factor of 2 in the implementation (but not in the plot), as also seen in the total force being a sum of the forces for F tot .Note that r is considered in the radial direction, effectively making this a 1D case.In (b) we have plotted F el + F ion for a range of particle radii to show their different y-intercept (where F el + F ion = 0).

Figure 3 .
Figure 3. Simulation run with 6000 particles of sizes 3.35μm in black and 3.65μm in red, ran for 1 second with radial geometry.(a) shows the initial positions of the particles, with (b) showing the final positions of the particles.(c) and (d) are the respective radial distribution of the species with b ¯being the mean radius for each species.

λ
De is the electron Debye length such that k

3. 1 . 1 . 1 ¯¯for
Discrete distribution For discrete binary distributions, the phase separation can simply be quantified by the mean radial position of each species b i ¯and computing the ratio b b 2 the respective species.In figure 3 we present final positions and radial distributions of simulated particles.In figures 3(c) and (d) the dashed vertical line shows the mean radii b ¯for initial and final time.As expected, we observed that the larger red particles move towards the edge and the smaller black particles move towards the center.In figure 4(a) the time evolution of b b 2 1

Figure 4 . 1 ¯¯is
Figure 4. Simulation with 6000 particles of binary size with mean radius a m 3.5m = ¯for different size disparity ò where a a a 2 1 = - ( ) ¯.In (a) the evolution of the system with the ratio of the mean radius b b 2 1 ¯¯is shown over time.In (b), the average value of b b 2 1

Figure 6 .
Figure 6.The slope β of the fitted line where a r C i b = + ¯, where a i¯is the average particle radii within one of 20 segment of radial distance from the center of the chamber r i where r is adapted to contain every single particle such that r r r min , max

Figure 7 .
Figure 7.The coupling parameter Γ versus standard deviation of particle size σ(a) and σ(a) −k , in (a) and (b) respectively, for all the distributions.The dashed line at Γ = 168 represents the phase transition between solid and fluid.The individual distribution types appears to be affected differently, with a log s G µ ( ) ( ) for uniform and discrete binary in (a), and a log k s G µ -( ) ( ) for Gaussian and binary Gaussian in (b), with k ≈ 0.3.Note that the discrete binary and binary Gaussian have different scaling for Γ and σ.

Figure 8 .
Figure 8.Translational pair correlation function g(r) for 4 sections radially.The rows indicate distribution types discrete binary, binary Gaussian, uniform and Gaussian, respectively.The column indicates size discrepancy with increasing σ(a) from left to right.On the left there are several peaks indicating a solid phase and long range order.Note the slight differences in binary and binary Gaussian, particularly for the blue line, indicating that discrete and continuous binary distributions do not behave exactly the same.

Figure 9 .
Figure 9. Orientational pair correlation function G 6 (r) for 4 sections radially.The rows indicate distribution types discrete binary, binary Gaussian, uniform and Gaussian, respectively.The column indicates size discrepancy with increasing σ(a) from left to right.On the left there are several peaks indicating a solid phase and long range order.Note the differences in binary and binary Gaussian, indicating that discrete and continuous binary distributions do not behave exactly the same.

r d 2 =
p , where r is the shell's radial position and d is the distance to the closest neighbour in the shell.Then, we have for r = 11 mm and d = 0.2 mm that N = 350.Then, the difference between expected distance to closest neighbor is approximately given by r

Figure 11 .
Figure 11.In (a): The average value of L • Ψ 6 over all particles versus the standard deviation of particle size for all the distributions.In (b) We observe a power law L a k 6 s Y µ -( • ) ( ) where k = 5.61, indicating that larger σ(a) exponentially degrades structure.Here, k is found by fitting L • Ψ 6 = cσ(a) −1/k such that L l o gc a log log k 6 1 s Y = + ( • ) ( ) ( ( )).

4. 2 .
Macroscopic structure and problems with structure factors While the coupling parameter in figure 7 has 2 different scalings for the Γ − σ(a) space, and at least 3 different slopes, the parameter L • Ψ 6 can sufficiently be predicted with one scaling and one slope as seen in figure 11(b).As such, we have macroscopic homogenity in the L • Ψ 6 − σ(a) across all distributions.However, in the Γ − σ(a) space we have strong heterogeneity with topological structures that depend on the individual particle size distrbutions.That being said, Γ can be predicted by σ(a), however that requires including the topological

Figure 12 .
Figure 12.The radial inter-particle repulsive force F repc , ion drag force F ion and electric field F el per particle.The dashed line is the ion drag force and dotted line is the electric field for particle radius a = 3.5μm.The discrete binary data is used in (a) and (b) while the binary Gaussian data is used in (c) and (d) both for the largest size disparity.

Figure 13 .
Figure 13.The distorted hexagonal lattice due to radial symmetry.Note the blue arrows which fails to intercept the third circle, and how the radial distance between shells shrinks closer towards the center.This helps visualize how the pair order parameters fail over long distances even though the hexagonal lattice is as good as it can be in radial geometry.

1 .
Phase separation, orientational and translational pair correlation functions can be predicted by σ(a) and are almost entirely independent of distribution type.All distribitions in the L • Ψ 6 − σ(a) space scale according to the same k = 5.61.

3 . 4 .
The coupling parameter space Γ − σ(a) scales differently depending on the distribition type.Here, Different size distributions have different final temperature, which in turn affects the coupling parameter.

Table 1 .
The set of selected radii for the dust particles for the different distributions in the 2D runs for Discrete Binary, Uniform, Gaussian and Binary Gaussian with the standard deviation σ as a fraction of the mean value μ a for the Gaussian distributions.The distributions are centered on a mean μ a such that μ a = 3.5 μm where Δa[μm] refers to the difference between the two means.The variations in column (a) and (d) are typically referred to as 'min' and 'max' respectively to indicate the degree of the size variation in figures.
a Here Δa refers to the difference between the largest and smallest particle.
[35,36]parameter[35,36]P H , given by screened away from the dust cloud such that n e is effectively zero and Z d = n i /n d .But, this assumes an infinitely large » ( )Where n i,0 is the undisturbed ion density i.e. the plasma as if there were no dust particles.Enforcing charge neutrality we have n i,0 = n e,0 + Z d n d , and we obtain P H ≈ 0.003 from the parameters in table 2, where n -• we obtain P H > > 1, resulting in electrons being

Table 3 .
Coupling parameter Γ for increasing size disparity from (a) to (d) respectively for the different size distributions.

Table 4 .
The size disparity ò and σ(a) for the discrete and Gaussian binary distributions.The distributions with 2 Gaussian peaks have a larger σ(a) despite the peaks' mean radii being equivalent to the discrete binary distribution.