Self-sustained non-equilibrium co-existence of fluid and solid states in a strongly coupled complex plasma system

A complex (dusty) plasma system is well known as a paradigmatic model for studying the kinetics of solid-liquid phase transitions in inactive condensed matter. At the same time, under certain conditions a complex plasma system can also display characteristics of an active medium with the micron-sized particles converting energy of the ambient environment into motility and thereby becoming active. We present a detailed analysis of the experimental complex plasmas system that shows evidence of a non-equilibrium stationary coexistence between a cold crystalline and a hot fluid state in the structure due to the conversion of plasma energy into the motion energy of microparticles in the central region of the system. The plasma mediated non-reciprocal interaction between the dust particles is the underlying mechanism for the enormous heating of the central subsystem, and it acts as a micro-scale energy source that keeps the central subsystem in the molten state. Accurate multiscale simulations of the system based on combined molecular dynamics and particle-in-cell approaches show that strong structural nonuniformity of the system under the action of electostatic trap makes development of instabilities a local process. We present both experimental tests conducted with a complex plasmas system in a DC glow discharge plasma and a detailed theoretical analysis.

A complex (dusty) plasma system is well known as a paradigmatic model for studying the kinetics of solid-liquid phase transitions in inactive condensed matter. At the same time, under certain conditions a complex plasma system can also display characteristics of an active medium with the micron-sized particles converting energy of the ambient environment into motility and thereby becoming active. We present a detailed analysis of the experimental complex plasmas system that shows evidence of a non-equilibrium stationary coexistence between a cold crystalline and a hot fluid state in the structure due to the conversion of plasma energy into the motion energy of microparticles in the central region of the system. The plasma mediated non-reciprocal interaction between the dust particles is the underlying mechanism for the enormous heating of the central subsystem, and it acts as a micro-scale energy source that keeps the central subsystem in the molten state. Accurate multiscale simulations of the system based on combined molecular dynamics and particle-in-cell approaches show that strong structural nonuniformity of the system under the action of electostatic trap makes development of instabilities a local process. We present both experimental tests conducted with a complex plasmas system in a DC glow discharge plasma and a detailed theoretical analysis.
Interest in model systems that permit a study of phase transitions in condensed matter started in the 1930s with the theory of the Wigner crystal 1 . Model systems provide significant insights for furthering our understanding of mechanisms of formation and growth of new phases and scenarios of phase transitions. In the scope of this problem, macroscopic systems with a high degree of non-ideality, such as colloidal suspensions in aqueous solutions 2 and dusty plasmas 3 , are of a special interest. Their properties can be studied in detail on a particle-resolved level.
In this paper, the focus is on dusty plasmas which are actively employed as a paradigm for studies of phase transitions 4,5 , collective wave excitations 6 , self-organization 7,8 and transport processes 9 . Micron-sized particles (dust) injected into a gas discharge plasma acquire a high negative electric charge of the order of 10 4 e by collecting the highly mobile electrons in the plasma. These highly charged dust particles can constitute a strongly coupled (non-ideal) component of the system and can form organized structures like crystals. These structures can be studied at the particle (kinetic) level by video-microscopic methods. Although dusty plasmas are often considered as toy models for the study of inactive condensed matter, under certain conditions it mimics active matter 10 and reveals similarity with such objects as cells in tissues 11 , suspensions of bacteria 12 and even the flockings of birds 13 . The common principle uniting all of these seemingly different systems is the specific mechanism of energy exchange between the system and the environment. Energy is injected into such systems at the level of each individual particle and is converted into the motion of particles. This leads to an intrinsic out-of-equilibrium behavior of the system and its departure from the action-reaction symmetry 14 at the microscopic level. The source of injected energy might be, for example, chemical interaction, as in self-phoretic Janus colloids 15 .
In dusty plasmas, there are several mechanisms that might make dust particles behave similar to agents in active matter. These mechanisms include the effects of ion and neutral shadowing 16 and the "rocket effect" 17  www.nature.com/scientificreports/ There are also experimental studies of Janus particles in a gas discharge plasma which become active due to the action of photophoretic force from the illuminating laser 18,19 . Another mechanism leading to energy conversion under certain conditions is the plasma-specific wake effect [20][21][22] . The wake effect arises due to the action of electric field that compensates for gravity acting on dust microparticles. This electric field also exerts a strong vertical ion flow which is perturbed by a highly charged particle of dust. The radial distribution of the electrostatic potential in the wake region around each dust particle is no longer spherically symmetric [23][24][25][26][27][28] . Asymmetry of the electrostatic potential and momentum exchange with the plasma flow lead to the nonreciprocal character of particle interactions in the dust subsystem. Particles absorb energy and momentum from the flowing plasma and then transport it from the local scale to the larger scales.
There are several examples of specific mechanisms in complex plasmas that might lead to an ordered-disordered structure coexistence in the dust subsystem although such coexistence has not been studied in detail before. Nosenko et al. 29 conducted an experiment where spontaneous formation of spinning particle pairs (torsions) was studied in a mono-layer complex plasma crystal by reducing the discharge power at several different values of a neutral gas pressure. The number of torsions is found to increase with a decrease in the discharge power. At lower gas pressures, the formation of torsions is preceded by a mode coupling instability (MCI). Development of MCI led to the observation of a crystal-fluid coexistence state in the monolayer. However, this state was not studied in detail: structural or thermodynamic parameter measurements were not carried out to delineate the transient nature of coexistence. In another experiment, Uchida et al. 30 produced a two-dimensional dust vortex flows in an ordered complex plasma structure. The vortex structures are produced by the effect of an asymmetry of ion drag force. The ordered structure is found to be disturbed by the dust flow, and disordered structure is produced in the vortex region which leads to an ordered-disordered structure coexistence in the system. Moreover, in a very recent molecular simulations study by Qiu et al. 31 , a two-dimensional Yukawa solid and liquid separation after the shock propagation is reported. After the propagation of compressional shocks, the structure and dynamics of the post-shock region are investigated. When the compressional speed is significantly higher than a first threshold value, the post-shock region melts completely. However, when this compressional speed is lower than a second threshold value, which is smaller than the first threshold, the post-shock region is found to be in the solid state. Whereas, at a compressional speed value between first threshold and second threshold, the post-shock region clearly exhibits the coexistence of the solid close to the compressional boundary and the liquid in the other part. Further, the averaged kinetic temperature in the post-shock region shows a spatial variation and is attributed to the dynamical heterogeneity of the 2D Yukawa systems.
In the present work we focus on the possible coexistence regime which is due to the plasma-specific mechanism of non-equilibrium melting in complex plasmas. This mechanism is caused by the wake effect and may arise mainly due to two types of instabilities. In multilayer systems, with at least two layers, the so-called "Schweigert" instability 32-34 develops in the form of driven horizontal oscillations of particles. It is due to the non-reciprocal interaction of dust particles located at different levitation heights. In a single-layered (mono-layer) quasi-2D system of dust particles another type of instability can arise due to the coupling of the horizontal motion of dust particles to their vertical motion. In this case of a mode coupling instability 35,36 particles oscillate with approximately equal amplitudes and frequencies in both the horizontal and vertical directions 27 . The conditions for the onset of MCI are determined by the strength of the vertical confinement, the density of the system and the neutral gas damping rate.
In most experimental works considering non-equilibrium melting of dust structures under the action of instabilities only melting of the entire structure is reported. The intermediate state where coexisting phases can be observed has not been studied experimentally. However, as soon as the threshold for the instability onset in dusty plasmas is density dependent in most cases, phase coexistence can be expected if the collection of particles is spatially nonuniform 37 . Spatial nonuniformity of structural and dynamical characteristics reveals in conditions of many experiments with dusty plasmas due to the action of central electrostatic confinement which keeps the system stable [38][39][40] . Such central trap can only be compensated by an inter-particle interaction force gradient and subsequent density gradient. The recent theoretical study of MCI-induced melting in the finite-size dust monolayer demonstrates that, if density gradient in the system is strong enough, then the peripheral part of the structure might stay ordered and coexist with the molten central region where MCI is active 41 . Moreover, Schweigert instability may also trigger a coexisting phase if the instability is driven locally in the system. This scenario demonstrates the applicability of phase coexistence concept to the collections of "active Brownian particles" in dusty plasmas where the particles in the locally molten region acquire energy from plasma through a non-equilibrium process. In the present paper, we present such a system where a non-equilibrium solid-fluid coexistence exists and the dust particles in fluid phase are driven active by the ion wake mediated non-reciprocal interaction.
We provide the detailed experimental and theoretical study of coexistence of solid and fluid states in a strongly coupled system of dust particles in a gas discharge plasma. In the experiment, the system of dust particles has complex geometry: multilayer 3D structure in the central region and single-layered quasi-2D structure at the periphery. Detailed profiles of kinetic temperature are given with particle-resolved precision. It is shown that a strong temperature gradient is present in the system: "hot" central 3D-liquid part coexists with "cold" peripheral 2D-crystalline part. The co-existing parts are in the state of dynamical equilibrium, thus the co-existence is stationary and self-sustained. The heat flux from the central part to the peripheral region is supported by the continuous energy input from the wake-mediated particle interaction. Theoretical description of the system is based on the multiscale approach to simulations: while the dynamics of dust particles is described using the Newton equation, the pairwise inter-particle interaction is calculated by the particle-in-cell method accounting ion-neutral collisions. In simulations the onset of instability in the central region along with the absence of instability at the periphery is demonstrated.  42 , which has an L-shaped vacuum vessel made of pyrex-glass consisting of two chambers as shown in Fig. 1a. The plasma and the complex (dusty) plasma are produced in the primary chamber, the auxiliary chamber provides connections to different subsystems for evacuating the chamber and producing and characterizing the plasma and the complex plasma. A base pressure of 0.1 Pa is initially achieved in the experimental chamber by a rotary pump, and Argon gas is then flushed several times through it to remove the impurities. The working gas pressure is finally set at ∼ 8 Pa using a mass flow controller. In the DPEx-II device, an asymmetrical electrode configuration consisting of a circular-shaped stainless steel anode, and a long tray-shaped cathode are employed. Fig. 1b represents the schematic diagram of this asymmetrical electrode configuration. A DC glow discharge Argon plasma is produced between the electrodes by applying a DC voltage of 450 V. The plasma current is estimated by measuring the voltage drop across a current limiting resistor of 2 k that is connected in series with the power supply. The plasma parameters such as the electron temperature, plasma density, plasma potential, etc., are measured using single and double Langmuir probes as well as emissive probes. For this discharge condition, the plasma density and the electron temperature are found to vary over a range of 0.8-2×10 15 m −3 and 2-4 eV, respectively. Mono-dispersive melamine-formaldehyde (MF) particles of diameter 10.66 μm are then introduced in the plasma by a dust-dispenser to form the complex plasma. These dust particles in the plasma environment get negatively charged by accumulating more electrons (due to their higher mobility) than ions and levitate in the cathode sheath region due to a balance between the downward gravitational force and the upward electrostatic force resulting from the cathode sheath. The dust particles get confined horizontally due to the sheath-electric field around a confinement ring placed on the cathode and form a circular-shaped mono-layer crystal . The micron-sized particles are illuminated by a green laser and a CCD camera is used to capture the Mie-scattered light. In our experiments, the width of the laser used for illuminating the particles is comparatively smaller than the inter-layer separation, which allows us to scan the particles layer by layer. However, the particles which move randomly due to their stochastic thermal motion in the liquid state are also considered in the analysis. A sequence of images is stored in a computer for further analysis. IDL and MATLAB-based software are employed for tracking the individual particles over time to study their dynamics. A hexagonally ordered complex plasma structure is observed when the neutral gas pressure and the discharge voltage are set at 8 Pa and 450 V, respectively. Production of such a stable mono-layer complex plasma crystal in the cathode sheath region of a DC glow discharge plasma was recently reported by Hariprasad et al. 39 . At that given discharge voltage, the plasma density decreases with a reduction of neutral gas pressure, which causes the sheath around the confining ring to thicken. As a result, the confinement area for the dust particles shrinks, which causes a few particles at the center of the mono-layer to move below the mono-layer as shown in Fig. 2b 38 . These particles then interact strongly with the particles that reside in the mono-layer and form a fluid-like structure in the central region. Thus, with the decrease of neutral gas pressure, a purely crystalline state changes to an ordered-disordered coexisting state, which remains stable over time. Fig. 2a shows a typical image of the coexistence of ordered-disordered states at a neutral gas pressure of ∼ 7.3 Pa. At that specific discharge condition, the central region of the complex plasma exhibits a disordered fluid-like structure surrounded by an ordered crystalline structure. It is also found that the particles in the center move randomly and don't possess an equilibrium position whereas they arrange themselves into an ordered state at the periphery and oscillate around their equilibrium position. The structural and the thermodynamical properties of these distinct regions of coexisting states differ significantly, and they will be delineated in detail in the upcoming sections.
Coupling parameter estimation. The Coulomb coupling parameter ( Ŵ ) is defined as the ratio of the potential energy of the dust particles with the kinetic energy 43,44 and it determines the phase state of the complex plasma system. In the past, many studies have been carried out to investigate the dependency of complex plasma system on the Coulomb coupling parameter 39,43,45 . Ikezi et al. 43 predicted that the Coulomb system exhibits an ordered structure when the Coulomb coupling parameter crosses a threshold value of Ŵ ∼ 172 if the definition   44 predicted a threshold value of the effective coupling parameter ∼ 106 using the definition via the mean inter-particle distance. The difference between the numerical values 172 and 106 for an unscreened system is explained by the ratio of the Wigner-Seitz radius to the mean inter-particle distance. Knapek et al. 45 introduced a new measurement technique to estimate the Coulomb coupling parameter and the dust temperature from the position information of dust particles. In this approach, dust particles are assumed to be distinguishable classical particles and obey Maxwell-Boltzmann distribution. In local equilibrium, when the interaction of these dust particles with the plasma as well as with individual neutral atoms is, on average, balanced by neutral friction, the dynamics of individual particles in the lattice can in principle be described by a Langevin equation 45 . In such a case, one can write down the probability distribution 45 for each lattice cell as, with T being the particle temperature, E the Einstein frequency and m the dust particle mass. The standard deviation of the velocity distribution and of the displacement distribution independently yield the dust temperature and the coupling parameter, respectively. The standard deviation of the velocity distribution is given by σ v = T m and the standard deviation of the displacement distribution is given by In the present set of experiments, the Coulomb coupling constant is estimated from the knowledge of inter-particle distance and standard deviation of displacement distribution as Ŵ eff = � σ r 2 .
In our experiments, we estimated the coupling parameter and the dust temperature by employing the above mentioned technique 39 . From the position information of the dust particles, the displacement distribution is obtained along with the standard deviation in both ordered and disordered phases 39 . Inter particle distances are calculated from the particle positions to estimate the coupling parameter 39 . The value of the effective coupling parameter in the ordered peripheral region is estimated to be ∼ 210, while in the disordered structure at the center, the value turns out to be ∼ 35. Both the values and the statement of phase coexistence in the observed system are in agreement with the threshold values obtained by Vaulina and Khrapak 44 and also with our earlier experiments 38,39 . The crystal-fluid phase coexistence of the complex plasma system is further confirmed by the Voronoi diagram based structural analysis.
Voronoi diagram based structural analysis. The coexisting structure is further analysed using a Voronoi diagram 46 , which segregates the plane with n points into convex polygons such that every polygon contains exactly one generating point, and each point in a given polygon comes closer to its generating point compared to others. This diagram provides the basic unit cell of the structure under investigation with the information of coordinate number and the defects present in the system. Fig. 2c depicts the Voronoi diagram of the phase coexisting structure of Fig. 2a. It is seen from the Voronoi diagram (depicted in Fig. 2c) that the outer portion of the structure consists of hexagonally oriented unit cells and forms a crystalline state. In contrast, the central portion of the Voronoi diagram consists of completely disordered polygons that are filled with defects and therefore essentially depict a fluid state. The amount of disorder and the number of defects become maximum at the center, and they reduce as one moves to the periphery. Figure 2d presents a zoomed view of the crystal-fluid coexisting boundary region in which the top and the bottom portion of the complex plasma correspond to crystalline and fluid states, respectively. There does not exist a rigid boundary that separates the coexisting phases; instead, the defects reduce towards the periphery. The transformation from the disordered fluid to the crystalline Non-equilibrium nature of the coexistence state. An important concern regarding a phase coexistence structure is the nature of its thermodynamic state. Generally, in equilibrium systems, the temperatures of the two coexisting phases remain the same 47 whereas, in a non-equilibrium phase coexistence state, the phases can exhibit different temperatures 48 . To determine the thermodynamic nature of our experimental phase coexistence state, we have estimated the temperatures at different locations of the dusty plasma to cover the regions of the two phases. For this, the trajectories of individual particles at various locations are examined from a sequence of recorded images in order to estimate their average velocities and the dust temperature is estimated as explained in the previous section. By assuming the dust particles obey Maxwell-Boltzmanian distribution, the average temperature of the dust particles is estimated from the experimentally obtained velocity distribution.
The full width at half maximum ( T = mσ v 2 , where, σ v is the standard deviation of the velocity distribution and m is the mass of the dust particles) of the distribution function gives the average temperature of the particles 39,45 . Figure 3 displays the spatial temperature profile of the complex plasma which has two coexisting phases. The '0' location refers to the center of the complex plasma. Figure 3 shows that the dust temperature has a maximum value ( ∼ 30 eV) in the fluid phase (at center), which then reduces by almost ten folds ( ∼ 3 eV) in the crystalline phase (at the periphery). The melting temperature is depicted in the Fig. 3 as a dotted line, to distinguish the coexisting phases more clearly. This drastic difference of temperatures between the two coexisting phases essentially implies the non-equilibrium nature of the system. The temperature profile along the crystal axis looks like the shape of a Gaussian distribution. This is because the dust density is higher at the center where the ion wake driven instability heats up the system more efficiently. The dust density reduces when one goes away from the center and it results in the deduction of the strength of instability and as well the particle temperature. This gives rise to the nonuniformity in the temperature profile. Moreover, the heat produced at the central region diffuses towards the periphery. At the crystal-fluid interface, the heat exchange continues and the crystal temperature rises even though there are no particles beneath. Such a smooth transition of dust temperature was also observed in laser heated dusty plasma systems in the past 49 . Interestingly, in our experiments the coexisting structure is found to be stable over time and both the phases continue to remain at different temperatures even if the system is left for a longer time. Since the complex plasma system is dissipative in nature, there has to be an energy source that heats the central region and keeps the system in a self-sustained non-equilibrium phase coexistence state. One possibility is the existence of an instability mechanism, such as a Schweigert instability 50,51 , that can heat the central portion of the complex plasma system through the ion wake mediated non-reciprocal interaction between the dust particles in the upper and lower layers. This will be briefly discussed in the later part of this paper.
In the past, Nosenko et al. 29 also conducted an experiment with a two-dimensional dust crystal in the sheath region of the plasma. By reducing the RF discharge power at high gas pressures they managed to change the strength of horizontal and vertical confinement. Due to this change, a few particles left the crystal plane, got trapped at the bottom and formed particle pairs (torsions) exhibiting either a circular or a more complex motion in the ion wake region. Formation of torsions finally led to the disordering of a large area of the crystal. They attributed this to a non-reciprocal Schweigert effect. They concluded that once a torsion was formed it became a source of disturbance in the plasma crystal. Through inter-particle interactions the energy of this disturbance caused by the flowing ions was redistributed into the surrounding lattice. At a particular low pressure of 12 Pa, they observed a mode coupling instability preceding the formation of torsions. MCI led to another melting scenario. Regarding these different scenarios they concluded that vertical pair formation and MCI are two competing scenarios for the mono-layer crystal response to a weakening of the vertical confinement. The value of neutral gas pressure seemed to be a key factor in determining which process took place first. These two cases can be easily distinguished by the significant difference in particle trajectories. In our experiments, the triggering www.nature.com/scientificreports/ mechanism for the melting of the central region is the Schweigert effect that is akin to the Nosenko's first set of observations. In contrast to the Nosenko et al. experiment, we do not observe any MCI induced melting. Under our conditions the melting and coexistence occur immediately after the formation of a multi-layered structure in the central region while the initial single-layered system is ordered. Also, we do not observe the formation of particle torsions as of Nosenko et al. even at lower gas pressure of ∼ 7 Pa.

Multiscale simulations of phase coexistence in complex plasmas
In order to study the underlying mechanism of phase coexistence described in the experimental section, we perform a series of MD simulations for the system of 2500 dust particles similar to the one observed in the experiment. The description of interaction between dust particles in MD simulations is based on the explicitly calculated distribution of electrostatic potential around dust particles. The prodecure of potential calculation allows to account for the wake effects in the system [20][21][22] . Many important plasma parameters required for the potential calculation, such as the temperature and concentration of electrons, the ion flow velocity and the electric field gradient,can be only approximately estimated from conducted experimental measurements. For this reason, the spatial distibution of electrostatic potential is calculated using several sets of parameters which fall into the experimentally estimated range. The procedure allowing to calculate the wake potential self-consistently via the PIC approach 52-55 is given in Supplementary Materials. The typical view of the spatial distribution of the electrostatic potential around a dust particle under considered conditions is shown in Supplementary Fig. S1. The comparison of potential profiles at different values of chosen parameters is given in Supplementary Fig. S2. The obtained collection of results for different wake potentials allows to identify the processes leading to nonequilibrium separation of phases in the experimental system. For MD simulations, mass m and diameter D of dust particles in the model system are chosen equal to the experimental values. The value of charge of a dust particle Q corresponds to the value of grain charge in the particular wake potential calculation. The time step of molecular dynamics equals 10 −3 /ω pd , where ω pd = Q 2 /md 3 is the plasma-dust frequency and d is the distance between dust particles.
We include the additional central potential into simulations that confines the likely charged particles and keeps the system stable: where r i = (x i , y i , z i ) is the radius-vector of an i-th particle, α and β are the horizontal and the vertical trap parameters, respectively. At each value of the particle charge, the values of trap parameters are found to obtain the same system size and geometry as in the experiment. In all simulations, α is of the order of 10 −2 cgs units, β is of the order of 1.0 cgs units. It is important to note that the strength of vertical confinement β plays a crucial role in the development of instabilities in a dusty plasma system 27,28 . For this reason, the value of β is varied in each simulation to observe the effect on the system properties.
The overall equation of motion that is solved numerically for the i-th dust grain has the form where the last two terms are mathematically equivalent to the action of Langevin thermostat at a room temperature 56 . The term −mγṙ i describes the friction of dust particles in the viscous gas environment. Its effect can be significant because the value of the damping factor γ can control the development of instabilities in the system. Dependence of the damping factor on the neutral gas pressure can be calculated using the following approximate formula 7 : where v n is the thermal velocity of neutral gas atoms. In order to study how instabilities develop in the simulated dust structures, the value of γ is varied around the value that corresponds to the pressure of neutral gas in the experiment. Each simulation with unique values of Q, α , β and γ starts with randomly positioned particles which interact via the reciprocal Yukawa potential. Under the action of confinement and mutual Yukawa interaction with each other, likely charged particles move to their equilibrium positions. Their kinetic temperature at this moment is equal to the one of the applied Langevin thermostat. Then the interaction of particles is substituted with the chosen wake potential calculated for the charge Q and the ion flow velocity v fl . System dynamics comes to the steady state, and then principal data is calculated. We use cutoff radius for the interaction 10 times higher than the average inter-particle distance in the structure.
The typical view of the simulated system is given in Fig. 4. It has the same radial size as the experimental structure and a similar value of inter-layer spacing in the central part. Depending on the value of β , the simulated structure can contain a different number of layers in the central region. At high values of β , the entire system is one-layered. With the decrease of β , additional layers form in the center of the structure and the single-layered peripheral subsystem decreases in size. This is an expected behavior for a dust system which transforms from a pancake-like structure into a spherical one at β/α = 1 57,58 . In the range β/α = 50 ÷ 300 the central region comprises two or three layers which is the case of interest similar to the experimental system. For this reason, β is varied in this range. www.nature.com/scientificreports/ Varying β and γ effects dynamics of the system and allows to observe phase co-existence in several structures with different grain charges Q. The important feature of the observed phase co-existence in all of the structures is the presence of strong temperature gradient between the center and the periphery of a structure. Under temperature, we mean the parameter of the velocity distribution which has a Maxwellian profile. The temperature gradient is due to the development of a wake-induced instability in the central region after the formation of the second layer. Intensive energy release has a local character and reaches the periphery of the system in the form of the heat flux. As the heat flux is not enough for the melting of the peripheral part, we observe non-equilibrium phase coexistence of a 3D central "liquid" and a 2D peripheral "crystal" in the simulated system.
The radial profiles of in-plane kinetic temperature of dust particles in the structure are given in Fig. 5 for the values Q = 15, 000 e and α = 0.035 cgs units. The temperatures are averaged over multiple runs. It can be seen that the increase of β leads to the decrease of temperature profile in the given system. The increase of β compresses the system in the vertical direction and strengthens the vertical confinement. The variation of γ effects the phase co-existence process as expected: the lower γ , the more intensive heating of the central region. With the decrease of γ , the phase boundary moves into the crystalline section of the structure until it melts entirely.
In order to demonstrate phase co-existence quantitatively, we rely on the parameter of inter-particle distance fluctuation (IDF) whose applicability is recently demonstrated for a dust monolayer 41 . This parameter can be used to define local phase state in a system where phases co-exist. Formulated for finite-size systems, IDF can be applied both to dust structures 59,60 and nanoparticles 61 . It has advantages over the popular Lindemann and local order parameters 62 in the considered case because it allows calculation for a system including both 3D and quasi-2D sections. At the same time, its application is simple: the subsystem is supposed to lose order when the numerical value of IDF is above 0.1.
As soon as our system contains 2500 particles, IDF can be used directly without modifications 5 . It is calculated by the formula: www.nature.com/scientificreports/ where r ij is distance between particles i and j; N is the number of particles used to calculate the parameter. Note that IDF is calculated locally, not over the entire structure but over small subsystems containing up to 15 particles. This allows to track phase state locally both in the central region and at the periphery. The radial profiles of IDF are shown in Fig. 6 both for the states with and without phase coexistence in the system with Q = 15, 000 e, α = 0.035 cgs units and β/α = 230 . The IDF parameter (called "Berry parameter" or even "Lindemann parameter" in a few works) was originally employed for nanoparticles of inert gas atoms. It was shown that when its value calculated for the entire particle reaches 0.1, it indicates melting of the structure. Later, the parameter was adapted for small clusters of Yukawa particles in the modified form via the variance of its block averaged value. Recently it was shown that the original form of IDF can be employed in the local approximation for a monolayer of dust particles. In this case, IDF is calculated for each separate crystalline cell and then averaged over the subsystem of interest. This local form of IDF has a nonuniform radial profile in a monolayer and roughly indicates melting of its shells at the same critical value of 0.1. In this work, this value is also applicable to the central section of the observed experimental system. The radial distribution function constructed separately for the central section of the structure indicates the loss of long-range order when the local value of IDF is above 0.1. The main advantage of IDF over the classical Lindemann parameter in the considered case is that it is calculated from relative particle positions and is not susceptible to fluctuations arising in the finite system. The type of instability that drives the observed phase co-existence is discussed in the next section.

Analysis of instabilities observed in simulations.
While geometry of the considered dusty plasma structures in MD simulations is identical, we observe different behavior of the simulated systems with different values of grain charge and wake potentials. The behavior depends on several factors. First, the value of grain charge determines the strength of interaction and the range of melting temperatures in the system. Second, the important factor is the principal nonhomogeneity of the structure: it is dense and two-layered in the central region and rarefied single-layered at the periphery. For this reason, necessary conditions for the instabilities that provide the energy input into the system are satisfied only in the bilayer central subsystem and have a local character. Finally, the type of instability that develops in the central region is determined by the inter-layer distance and the structure of a wake potential.
The stability principles for bilayer complex plasmas are well studied in case of extended systems 27,28 . There are two types of instabilities that might arise in the bilayer system. The dynamical instability, also known as the Schweigert instability 32-34 , is of the oscillatory type. It arises due to a wave mode that is pumped by the energy of the ion flow and grows exponentially until damping limits its growth. This mode can be suppressed by sufficiently strong neutral gas friction. It has an important fingerprint that allows to identify its development: in case of the Schweigert instability, particles mainly oscillate in the in-layer direction. Vertical oscillations become negligibly small compared to the horizontal ones. Oscillations in the horizontal direction occur at a frequency which typically lies in the range 0.6ω pd − 1.1ω pd for the conditions of dusty plasma experiments.
Another type of instability that might arise in the bilayer system is the structural instability which is not of the oscillatory type 28 . It is due to the state of neutral equilibrium when the equilibrium positions of particles no longer correspond to the ground state of the system. The extended motion of particles is present in the structure for any damping. Note that both the dynamical and the structural instability can be mixed.
In our MD simulations, we observe that both types of instabilities can develop locally only in the bilayer central region of the structure. To identify the type of instability, we employ the following technique. Starting from the phase co-existing state, the friction coefficient γ is increased to cool down the central region of the system. In this case we find that the disordered, flow-like motion of particles is present in the bilayer central region even when the friction coefficient γ is many times higher than in the experiment. There is no fingerprint of an oscillatory instability type in the spectrum of particle oscillations.
In most simulations it is the Schweigert-type instability that is responsible for the heating of the central region. An example of the oscillation spectrum for the parameters Q = 15, 000 e , α = 0.035 cgs units, β/α = 150 , γ = 31.0 s −1 is given in Fig. 7. The main peak of in-plane oscillations and the ratio of horizontal to vertical amplitudes identifies the dynamical instability of Schweigert type. It develops locally in the multi-layered central region and fades behind the interface between the 3D and the 2D sections of the structure. The single-layered peripheral subsystem might only become unstable due to the mode coupling instability which does not develop in the considered conditions. For this reason, the following physical picture of the phase co-existence process can be formulated.
When particles in the central region organize into the second layer, under certain conditions the structural or the Schweigert-type instability arises. Each dust particle starts to act as an active agent and converts energy of the flowing plasmas into its own kinetic energy. The conversion of energy occurs only in the bilayer central region where the conditions for the development of instability are satisfied. Due to the conversion of energy, the central region melts. The peripheral subsystem where no instability is present is heated only by the heat flux from the central region. The flux of heat increases kinetic energy of particles at the periphery but it is not enough to melt the subsystem over a wide range of system parameters. In this case, we observe solid-liquid phase co-existence that is supported by active agents located in the center of the structure.

Comparison between theory and experiment and discussion
The employed simulation technique allows to obtain good agreement with the experimentally measured profile of kinetic temperature in the structure. The profile is demonstrated in Fig. 8 for the estimated values of system parameters Q = 15, 000 e , α = 0.035 cgs units, β/α = 230 , γ = 1.8 s −1 . This value of γ equals to the experimental value calculated by the formula (4) within the experimental error. The wake potential in this case is calculated for the value of ion flow velocity v fl = v B which lies in the experimentally estimated range. Note that we do not expect the complete self-consistency of the developed multiscale approach due to the approximations employed in the numerical model. In order to improve the results, one should calculate the charge of each dust particle selfconsistently with the dynamics of all dust particles. This is a resource-intensive task for the system of 2500 dust grains.
Inspired by the observation of phase co-existence in our experiment, it is tempting to compare some of the features of our complex plasma system to those of the active matter models. In our case, the principal possibility of phase co-existence arises from the nonuniformity and thermodynamic openness of the system. Due to its finite size and action of the confinement, it has a bilayer structure in the central region and a monolayer structure at the periphery. As our MD simulations show, in the central region conditions for the development of wakeinduced instabilities are satisfied. In presence of the instability, each charged dust particle in the central region can act as an active agent with the ion wake beneath such a dust particle serving as a source of input energy on a micro-scale. The interaction between the dust particle and the ion wake breaks the symmetry in the force www.nature.com/scientificreports/ balance, leading to a non-reciprocal interaction between the dust particles in the upper and lower layers. This non-reciprocal nature of interaction acts as a heat source, which helps to transfer the energy from the plasma species to the dust particles. The energy gained by the dust particles in the form of kinetic energy is dissipated by the background neutrals. Akin to an active mater, the central part of the phase co-existing complex plasma system gets a continuous supply of heat energy from the particle beneath, which keeps this region melted. The peripheral region of the co-existing state, on the other hand, does not have such an energy source due to the absence of arising instability. As a result, it always remains in the crystalline phase. The heat energy also gets transferred from the center to the periphery along the longitudinal direction. The neutrals surrounding the crystalline state act as a heat sink and take away the excess energy from the particles in the crystalline phase.
Assuming that the action of the heat source in the single-layered peripheral section can be neglected, the heat transport equation in this region has the form div(κ∇T) = 2γ nTk B , where κ is the thermal conductivity, n is the areal number density of dust particles. The experimentally measured temperature profile in the peripheral region is in a good agreement with the numerical solution of this equation. We estimate the value of the thermal diffusivity χ = κ/cn , where c is the specific heat of the dust subsystem. The estimated value is χ = 5 mm/s 2 which is close to the reported value 9 mm/s 2 in a 2D dusty plasmas near melting 49 .
Our present experimental observations could serve as a paradigmatic example of a complex plasma system exhibiting dynamical features of an active matter with a stationary co-existence of solid and liquid phases. The proposed multiscale numerical model allows to simulate the process of phase co-existence in the system of active dust particles starting from the plasma microscale up to the macroscale where the dust system evolves.

Conclusion
We experimentally observed a paradigmatic example of stable non-equilibrium coexistence of solid and fluid states in a complex plasmas system. The non-equilibrium nature of observed coexistence is confirmed by the measured nonuniform spatial profile of temperature in the structure. We conducted multiscale simulations of the system based on the kinetic equation for the calculation of electrostatic potential and on the Newton equation for the dynamics of dust particles. It is demonstrated that the coexistence is driven by the wake effect leading to arising of instabilities in the system. Due to structural nonuniformity under the action of the electrostatic trap, instabilities in the system develop locally. Our results prove the possibility of phase coexistence scenario in onecomponent plasmas which is often considered as a model system for classical condensed matter. We observe similarity between systems of active particles and dusty plasmas under certain conditions. The demonstrated coexistence of ordered and disordered states can be considered as the coexistence of a non-activated subsystem which is in thermal equilibrium with the neutral gas and a subsystem activated by per-particle energy exchange with surrounding plasma. We hope our experimental observations along with the numerical modelling results can trigger further investigations into this intriguing state of phase coexistence and the acquired 'active' nature of a dusty plasma system under certain conditions.

Data availability
The datasets used and/or analysed during the current study available from the corresponding author on reasonable request.