1. Introduction
One of the features of dusty plasma, which distinguishes it from the ordinary three component plasma, is that due to the heavy mass of the dust, there is a vast separation of time scales in the dynamics of the dust and that of the background plasma. This vast separation of time scales has an interesting consequence that the dust component can be confined in a state of thermodynamic equilibrium in a small volume within a sufficiently large plasma background, by external fields. This remarkable feature is shared by only one other system, i.e. non-neutral plasma which is confined in the state of thermodynamic equilibrium by external magnetic fields (Davidson Reference Davidson1990). In this paper, we construct a general thermodynamical model of dusty plasma where the dust is confined in a small volume within a sufficiently large plasma background, by external fields. The model is solved analytically in the mean field limit and various processes for the gaseous phase of dust, e.g. isothermal/adiabatic/constant internal energy expansion/compression, free expansion of dust, specific heat and dispersion of acoustic waves, etc. are studied.
While defining thermodynamic processes in dusty plasma, the role of plasma background must be carefully examined. The reason for this is that in these plasmas, the dust component is coupled with the background electrons and ions through quasi-neutral electric fields. Thus, when the dust component is compressed/expands or moves, for example in dust acoustic wave (DAW) or shocks, then the background plasma is also perturbed. The extra heat energy generated due to such perturbations of the background plasma must be properly accounted for while defining thermodynamic processes in dusty plasma. A somewhat truncated version of this model was described earlier to show the interconversion of plasma heat into work and vice versa via cold dust (Avinash Reference Avinash2010a,Reference Avinashb; Avinash & Kaw Reference Avinash and Kaw2014). However, since the dust temperature 
${T_d}$ was taken to be zero in the model, the thermodynamic processes involving the dust internal energy, pressure, enthalpy and other related thermodynamic parameters could not be described properly. In this paper, the thermodynamic model with finite dust temperature is discussed.
Hamaguchi & Farouki (Reference Hamaguchi and Farouki1994) and Farouki & Hamaguchi (Reference Farouki and Hamaguchi1994) in their seminal works, have also proposed a thermodynamic model of dusty plasma to calculate dust correlation effects and the solid–liquid melting boundary. In this model, the dust and the plasma background occupy the same volume. The negative charge of the dust is confined or neutralized by the cohesive plasma background. This, however, is not realistic. In dusty plasma experiments the dust is confined in a small volume within the plasma with the help of external electrostatic (ES) fields (Barkan & Merlino Reference Barkan and Merlino1995; Trottenberg, Block & Piel Reference Trottenberg, Block and Piel2006; Pilch et al. Reference Pilch, Piel, Trottenberg and Koepke2007; Thomas Reference Thomas2010). Our model takes into account this confinement of dust within the background plasma. There are some other important differences in the results of our model and those of the Hamaguchi–Farouki (HF) model which will be discussed later in the paper.
2. General formulation of the model
Our model consists of an ensemble of 
${N_d}$ identical point dust particles carrying a constant negative charge 
${q_d}$, having pressure and temperature 
${P_d},{T_d}$, respectively. The dust cloud is assumed to be confined in a small volume 
${V_d}$ by external fields 
${P_{\textrm{ext}}}({P_d} = {P_{\textrm{ext}}})$, within a weakly coupled, statistically averaged, plasma background, having 
${N_e}$ and 
${N_i}$ numbers of electrons and ions, in volume V where 
$V \gg {V_d}$. Typically, in experiments 
${V_d}/V \approx {10^{ - 3}}{-}{10^{ - 5}}$ (Trottenberg et al. Reference Trottenberg, Block and Piel2006; Fisher et al. Reference Fisher, Avinash, Thomas, Merlino and Gupta2013) which justifies the assumption of small 
${V_d}$ in the present model.
The electron and ion temperatures are denoted by 
${T_e}$ and 
${T_i}$, respectively. In thermal equilibrium the electron and the ion densities are given by Boltzmann's relations 
${n_i} = {n_0}{\rm exp} ( - q\varphi /{T_i}),\;{n_e} = {n_0}{\rm exp} (q\varphi /{T_e})$ where φ is the electrostatic potential within the dust cloud. Sufficiently away from the cloud φ → 0 and 
${n_i} = {n_e} = {n_0}$. It should be noted that this assumption is asymptotically valid in the limit 
${V_d}/V \to 0$, i.e. in the limit of dust cloud embedded in an asymptotically infinite plasma background. Another point to be noted is that the volume V and 
${V_d}$ can be varied independently in the present model, which is consistent with dusty plasma experiments (Barkan & Merlino Reference Barkan and Merlino1995; Trottenberg et al. Reference Trottenberg, Block and Piel2006; Pilch et al. Reference Pilch, Piel, Trottenberg and Koepke2007; Thomas Reference Thomas2010). This is different from the HF model where the dust density is assumed to be proportional to the average plasma density (Farouki & Hamaguchi Reference Farouki and Hamaguchi1994). In the last section, and further in appendix B, we will discuss these issues and the relationship of the present model and the HF model in detail.
The dust and the plasma are usually immersed in the background gas filling the plasma chamber. The collisions of the electrons, the ions and the dust with neutral atoms of the back ground gas help to regulate the temperatures of the dust and the plasma. Thus, the neutral background acts as the heat bath (Quinn & Goree Reference Quinn and Goree2000a). In dusty plasmas, the dust is heated by a number of processes, e.g. Brownian motion due to the neutral gas, fluctuations due to the electric field or the dust charge and it is cooled mainly due to the drag by neutrals (Quinn & Goree Reference Quinn and Goree2000a). The typical dust neutral collision frequency 
${\nu _{\textrm{dn}}}$ (
${\propto} {P_g}$ neutral gas pressure) in dusty plasma experiments ranges from a few Hz at low gas pressures to ≥100 Hz at high pressures (Quinn & Goree Reference Quinn and Goree2000b). Thus, if 
${\tau _p} \gt {\tau _{\textrm{dn}}}$ (
${\tau _p}$ is the time scale of the thermodynamic process and 
${\tau _{\textrm{dn}}} = 1/{\nu _{\textrm{dn}}}$ is the time scale of the dust neutral collision), which would typically be true in experiments with high neutral pressure (100–200 mTorr), then in such cases, the dust is in good thermal contact with the neutral bath. In this regime, the isothermal approximation for dust with 
${T_d} \approx {T_n}$ (
${T_n}$ is the temperature of neutrals) is appropriate (Thomas Reference Thomas2010). On the other hand, in experiments with low neutral pressures and relatively rare collisions with neutrals (≤10 mTorr) 
${\tau _p}$ typically is less than 
${\tau _{\textrm{dn}}}$. In such cases, the thermal contact of dust with the neutral bath is rather weak (Pilch et al. Reference Pilch, Piel, Trottenberg and Koepke2007). In this regime, the dust may be treated as isolated and an adiabatic approximation with variable temperature for the dust is appropriate. The thermal conductivity of electrons and ions in the background plasma is always very large. The ion heat is lost mainly through collisions with neutrals or due to the diffusive losses to the walls of the plasma vessel on time scale of approximately 1 to 10 
${\rm \mu}$s. The electron heat is lost on an even shorter time scale. Thus, for the plasma background, on the time scale 
${\tau _p}$, we assume an isothermal approximation with fixed electron and ion temperature. Hence, apart from the neutral background, the dust is also coupled to the plasma background, with which it exchanges heat via electric fields.
We start our calculations by expressing the energy conservation for quasi-static (
${\tau _p} \gt {\tau _{\textrm{relax}}}$ where 
${\tau _{\textrm{relax}}}$ is relaxation time scale) work done by 
${P_{\textrm{ext}}}$ on 
${V_d}$ given by
\begin{equation}\Delta {Q_e} + \Delta {Q_i} + \Delta {Q_d} = \Delta U + {P_d}\Delta {V_d}.\end{equation}
In this equation 
$\Delta {Q_d} = {T_d}\Delta {S_d},\Delta {Q_e} = T{}_e\Delta {S_e},\Delta {Q_i} = {T_i}\Delta {S_i}$ where 
$\Delta {Q_d}$ is the heat exchanged with dust heat bath, while 
$\Delta {Q_i},\Delta {Q_e}$ denote the heat exchanged with the plasma background and V is held constant. The entropies of the dust, electrons and ions are given by 
${S_d},{S_e}$ and 
${S_i}$, respectively, and U is the internal energy of the composite system of the dust and particles. The expressions for plasma entropies and U are given by (Hamaguchi & Farouki Reference Hamaguchi and Farouki1994; Avinash Reference Avinash2010a)
\begin{gather}U = \frac{3}{2}({N_e}{T_e} + {N_i}{T_i} + {N_d}{T_d}) + \frac{1}{2}\int {\rho \varphi \,\textrm{d}V} - \frac{{q_d^2}}{{8{\rm \pi} {\varepsilon _0}}}\sum\limits_{j = 1}^{{N_d}} {\int {\frac{{\delta (\boldsymbol{r} - \boldsymbol{r}_{\boldsymbol{j}})}}{{|\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{j}}}|}}} \,\textrm{d}V,}\end{gather}
\begin{gather}\left. {\begin{array}{c@{}} {{S_\alpha } = \tfrac{3}{2}{N_\alpha } - \int_V {{n_\alpha }[(\ln {n_\alpha }\Lambda_\alpha^3) - 1]} \,\textrm{d}V}\\ {{\Lambda_\alpha } = {{\left( {\dfrac{{{h^2}}}{{2{\rm \pi} {m_\alpha }{T_\alpha }}}} \right)}^{1/2}}} \end{array}} \right\},\end{gather}
where α denotes electrons or ion, and ρ and φ are the local charge density and the electrostatic potential, while h denotes the Planck constant. In the isothermal equilibrium of the background and the approximation 
$q\phi /{T_\alpha } \lt 1$ (justified later), the electron and ion densities are given by the linearized Boltzmann relations 
${n_i} = {n_0}(1 - q\varphi /{{T}_{i}}), {n_e} = {n_0}(1 + q\varphi /{{T}_{e}})$. The ES potential φ in the dust cloud can be obtained by solving the corresponding Poisson's equation given by
\begin{equation}{\varepsilon _0}{\nabla ^2}\varphi = \rho ={-} {q_s}\sum\limits_{j = 1}^{{N_d}} {\delta (\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{j}}}) - {\varepsilon _0}\varphi /\lambda _d^2} ,\end{equation}
where 
$1/\lambda _d^2 = ({q^2}{n_0}/{\varepsilon _0})(1/{T_e} + 1/{T_i})$, and we have substituted linearized Boltzmann relations for electrons and ion densities. The solution of (2.4) is given by
\begin{equation}\varphi ={-} \frac{{{q_d}}}{{8{\rm \pi} {\varepsilon _0}}}\sum\limits_j {\frac{{{\rm exp} ( - |\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{j}}}|/{\lambda _d})}}{{|\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{j}}}|}},}\end{equation}
where the index 
$j = 1 \ldots {N_d}$. As shown in appendix A (Hamaguchi & Farouki Reference Hamaguchi and Farouki1994; Avinash Reference Avinash2010a) using (2.2)–(2.5) we derive the following expressions for the internal energy and the heat exchanged with the isothermal plasma background:
\begin{gather}\begin{array}{ccccc} U & = \dfrac{3}{2}({N_d}{T_d} + {N_e}{T_e} + {N_i}{T_i}) - \dfrac{{q_d^2{N_d}{\kappa _d}}}{{8{\rm \pi} {\varepsilon _0}}} + \dfrac{{q_d^2}}{{8{\rm \pi} {\varepsilon _0}}}\sum\limits_i {\sum\limits_{j \ne i} {\dfrac{{{\rm exp} ( - |{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|/{\lambda _d})}}{{|{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|}}} } \\ & \quad - \dfrac{{q_d^2}}{{16{\rm \pi} \varepsilon {\lambda _d}}}\sum\limits_i {\sum\limits_{j \ne i} {{\rm exp} ( - |{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|/{\lambda _d})} } , \end{array}\end{gather}
\begin{gather}\begin{array}{ccccc} {T_i}{S_i} + {T_e}{S_e} & = \dfrac{3}{2}({N_e}{T_e} + {N_i}{T_i}) - \sum\limits_\alpha {{T_\alpha }{N_\alpha }[(\ln {n_0}\varLambda _\alpha ^3) - 1]} \\ & \quad - \dfrac{{q_d^2}}{{16{\rm \pi} \varepsilon {\lambda _d}}}\sum\limits_i {\sum\limits_{j \ne i} {{\rm exp} ( - |{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|/{\lambda _d})} } , \end{array}\end{gather}
where 
${\kappa _d} = 1/{\lambda _d}$ and indices 
$i,j = 1 \ldots {N_d}$. The second term on the right-hand side of (2.7) is due to the uniform isothermal plasma background. The energy conservation equation for the dust component can be constructed from (2.1) as
\begin{equation}\Delta {Q_d} = {T_d}\Delta {S_d} = \Delta {U_d} + {P_d}\Delta {V_d},\quad {U_d} = U - {T_e}{S_e} - {T_i}{S_i},\end{equation}
where 
${U_d}$ is the effective internal energy of the dust component. The expression for 
${U_d}$ can be obtained by eliminating 
$- ({T_i}{S_i} + {T_e}{S_e})$ and U from (2.6) and (2.7) to give
\begin{equation}\begin{array}{ccccc} {U_d} & = \dfrac{3}{2}({N_d}{T_d}_{}) + \dfrac{{q_d^2}}{{8{\rm \pi} {\varepsilon _0}}}\sum\limits_i {\sum\limits_{j \ne i} {\dfrac{{{\rm exp} ( - |{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|/{\lambda _d})}}{{|{\boldsymbol{r}_{\boldsymbol{i}}} - {{\boldsymbol{r}}_{\boldsymbol{j}}}|}}} } \\ & \quad + \sum\limits_\alpha {{T_\alpha }{N_\alpha }[(\ln {n_0}\varLambda _\alpha ^3) - 1]} - \dfrac{{q_d^2{N_d}{\kappa _d}}}{{8{\rm \pi} {\varepsilon _0}}}. \end{array}\end{equation}From (2.8a,b), the corresponding expressions for Helmholtz energy, pressure, entropy and enthalpy of the dust component can be calculated as
\begin{equation}F{}_d = ({U_d} - {T_d}{S_s}),\quad {P_d} ={-} {\left. {\frac{{\partial F}}{{\partial {V_d}}}} \right|_{{T_d}}},\quad {S_d} = {\left. { - \frac{{\partial {F_d}}}{{\partial {T_d}}}} \right|_{{V_d}}}, \quad H_{d} = (U_{d} + P_{d} \, V_{d}).\end{equation}
From these expressions the thermodynamic variables of the dust can be calculated as follows. The effective internal energy of the dust 
${U_d}$ can be obtained from (2.9). The free energy 
${F_d}$ can be directly obtained from 
${U_d}$ by integrating the thermodynamic relation 
${U_d}/T_d^2 ={-} \partial /\partial {T_d}({F_d}/{T_d})$. From 
${F_d}$, the pressure, the entropy and the enthalpy of the dust can be calculated from (2.10a–d). Typically, these quantities will have the usual thermal component and an excess 
${X_{\textrm{ES}}}$ due to electrostatic contributions. Since the last two terms in (2.9) depend on 
${T_e},{T_i}$ and V, which are constant, hence 
${U_d} = {U_d}({T_d},{V_d})$. Thus 
${U_d}$ and its 
${X_{\textrm{ES}}}$ can be calculated for given values of dust density and temperature from Molecular Dynamic (MD) simulations involving dust particles alone (Farouki & Hamaguchi Reference Farouki and Hamaguchi1994; Hamaguchi & Farouki Reference Hamaguchi and Farouki1994). In the next section, we show that in the thermodynamic limit, 
${X_{\textrm{ES}}}$ can be calculated analytically.
In the present model, the assumptions of constant dust charge 
${q_d}$ and 
$q\varphi /{T_e},q\varphi /{T_i} \lt 1$ used in the linearization of the Boltzmann response are valid in the limit 
${q_d}{n_d}/2q{n_0} \lt 1$ (Havnes et al. Reference Havnes, Goertz, Morfill, Grun and Ip1987; Goertz Reference Goertz1989), which is well-satisfied for typical dusty plasma experimental parameters. The dust temperature 
${T_d}$, on the other hand, is limited by the dust confinement. In dusty plasma experiments, the dust is electrostatically confined in a parabolic potential well 
${\varphi _{\textrm{Ext}}} \propto {r^2}$, which is nearly zero in the centre and rises monotonically to value 
${\varphi _E}$ at the edge of the confinement at 
$r = {r_E}$. Hence the condition for the confinement of dust in the ES potential well is 
${T_d} \lt {q_d}{\varphi _E}$.
3. Mean field limit
In the mean field limit, the dust is in the rare density gaseous phase. In this phase, the correlations are exponentially weak and the dust can be smeared out into a homogenous uncorrelated fluid. As will be shown shortly, in this limit there is a mean electrostatic field within the dust cloud. The mean field limit is given by 
${q_d} \to 0$, 
${N_d} \to \infty$, 
${q_d}{N_d} \ne 0$ and the double summation is replaced by smooth integrations as 
$\sum\nolimits_i {\sum\nolimits_j { = {n_d}{N_d}\int {\textrm{d}{V_d}} } }$ in (2.6) and (2.7). Carrying out double summations according to this prescription we obtain (Avinash Reference Avinash2010a)
\begin{equation}\begin{array}{ccccc} & \dfrac{{q_d^2}}{{8{\rm \pi} {\varepsilon _0}}}\sum\limits_i {\sum\limits_{j \ne i} {\dfrac{{{\rm exp} ( - {\kappa _d}|{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|)}}{{|{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|}}} = } \dfrac{{q_d^2}}{{16{\rm \pi} {\varepsilon _0}}}\sum\limits_i {\sum\limits_{j \ne i} {{\kappa _d}{\rm exp} ( - {\kappa _d}|{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|)} } \\ & \quad = \dfrac{{q_d^2N_d^2}}{{{q^2}n}}\dfrac{{{T_e}{T_i}}}{{{V_d}({T_e} + {T_i})}}. \end{array}\end{equation}Thus the dust internal energy, in the thermodynamic limit, is given by
\begin{equation}{U_d} = \frac{3}{2}{N_d}{T_d} + \frac{{q_d^2N_d^2{T_e}{T_i}}}{{{q^2}n{V_d}({T_e} + {T_i})}} + \sum\limits_\alpha {{T_\alpha }{N_\alpha }[(\ln {n_0}\varLambda _\alpha ^3) - 1]} ,\end{equation}
where 
$n = 2{n_0}$ is the pristine plasma density. In (3.2), the second term is due to the mean ES field in the dust cloud. Thus, in the thermodynamic limit, though the correlations are weak, dust particles still interact through the mean ES field. This is different from ideal gas where particles do not interact at all.
Substituting 
${U_d}$ in the thermodynamic relation 
${U_d}/T_d^2 ={-} \partial /\partial {T_d}({F_d}/{T_d})$ and integrating gives the Helmholtz free energy, pressure, entropy of the dust as
\begin{gather}{F_d} = {T_d}{N_d}[\ln ({n_d}\varLambda _d^3) - 1] + \frac{{q_d^2N_d^2{T_e}{T_i}}}{{{q^2}n{V_d}({T_e} + {T_i})}} + \sum\limits_\alpha {{T_\alpha }{N_\alpha }[(\ln {n_0}\varLambda _\alpha ^3) - 1],}\end{gather}
\begin{gather}{P_d} = \frac{{{N_d}{T_d}}}{{{V_d}}} + \frac{{q_d^2N_d^2{T_e}{T_i}}}{{{q^2}nV_d^2({T_e} + {T_i})}},\quad {S_d} = \tfrac{3}{2}{N_d} - {N_d}[\ln ({n_d}\varLambda _d^3) - 1].\end{gather}
In (3.4a), the second term gives the excess pressure 
${P_{\textrm{ES}}}$ due to mean electric fields in the confined dust cloud. This electrostatic pressure 
${P_{\textrm{ES}}}$ (Avinash Reference Avinash2010a,Reference Avinashb) can be dominant for typical experimental parameters.
The dust pressure 
${P_d}$ was experimentally measured by Fisher et al. (Reference Fisher, Avinash, Thomas, Merlino and Gupta2013). The value of 
${P_d}$ was found to be much greater than the thermal pressure 
${n_d}{T_d}$. In fact, the experimentally measured value of 
${P_d}$ was found to be within a factor of order unity of the value predicted by the ES part of total pressure 
${P_{\textrm{ES}}}$ given in (3.4a). The ES dust pressure 
${P_{\textrm{ES}}}$ was further confirmed in the experiments by Williams (Reference Williams2019) and MD simulations (Shukla et al. Reference Shukla, Avinash, Mukherjee and Ganesh2017). From (3.4b) it is seen that in the gaseous limit, the ES contributions to the dust entropy are zero.
4. Thermodynamic processes
In this section we define thermodynamic processes, specific heat, dispersion of acoustic waves and free expansion involving dust.
4.1. Isothermal process
In the case where the time scale of the dust process (e.g. expansion/compression), 
${\tau _P}$ is larger than the dust neutral collision time scale 
${\tau _{\textrm{dn}}}$, the dust is strongly coupled with neutrals. In this case, the dust temperature 
${T_d}$ is regulated by the neutral bath and can be taken to be constant 
${T_{d0}}({\approx} {T_n})$. This would typically be true in experiments with high gas pressure (Thomas Reference Thomas2010). Thus, in the high gas pressure regime, the equation of state is given by
\begin{equation}{P_d} = \frac{{{N_d}{T_{d0}}}}{{V_d^{}}} + \frac{{q_d^2N_d^2{T_e}{T_i}}}{{{q^2}n({T_e} + {T_i})V_d^2}}.\end{equation}
In this case, the thermal energy of the dust component 
$3{N_d}{T_{d0}}/2$ is constant. However, the internal energy 
${U_d}$ which contains ES contributions as well, and is defined in (3.2), is not constant. Processes where 
${U_d}$ is constant will be defined later. For quasi-static expansion/contraction of the dust volume 
$\Delta {V_d}$, the energy conservation can be re-expressed as 
${T_d}\Delta {S_d} + \Delta ({T_i}{S_i} + {T_e}{S_e}) = {P_d}\Delta {V_d}$, which shows that the heat taken from the plasma bath maintains an isothermal background with given electron and ion temperatures. The dust then expands/compresses at constant 
${T_d}$ in this isothermal background. The change in the dust entropy is given by 
$\Delta {S_d} = {n_d}\Delta {V_d}$ while the change in the plasma entropy (ES part), in the simple case where 
${T_e} = {T_i} = T$ is 
$\Delta S = q_d^2n_d^2\Delta {V_d}/2{q^2}n$. In terms of coupling parameters 
${\varGamma _0},{\kappa _0}$ which, for given 
${N_d}$, are related to the two constants 
${T_{d0}},{V_{do}}$ by the relations 
${\varGamma _0} = q_d^2/4{\rm \pi} {\varepsilon _0}{a_0}{T_{d0}},{\kappa _0} = {a_0}/{\lambda _d}$ where 
${a_0} = {(3{V_{d0}}/4{\rm \pi} {N_d})^{1/3}}$, the equation of state is given by
\begin{equation}{\bar{P}_d} = \frac{1}{{{{\bar{V}}_d}}} + \frac{3}{2}\frac{{{\varGamma _0}}}{{\kappa _0^2}}\frac{1}{{\bar{V}_d^2}},\end{equation}
where we have normalized pressure with 
${n_{d0}}{T_{d0}}$ and 
${V_d}$ by 
${V_{d0}}$. In the weak correlation regime the ratio of the ES potential energy to the average kinetic of dust is given by 
$\varGamma _0^\ast = {\varGamma _0}{\rm exp} ( - {\kappa _0}) \lt 1$. In figure 1 we show a dust isotherm for 
${\varGamma _0} = 1,{\kappa _0} = 1.5,\varGamma _0^\ast = 0.2$. The isothermal equation of state given in (4.2) has been verified in recent MD simulations (Shukla et al. Reference Shukla, Avinash, Mukherjee and Ganesh2017). In the large 
${V_d}$, limit, 
${\bar{P}_d} \propto 1/{\bar{V}_d}$ as an ideal gas, however, for smaller dust clouds where ES pressure dominates, 
${\bar{P}_d} \propto 1/\bar{V}_d^2$.
Figure 1. Isothermal process (solid), adiabatic (dash) and constant internal energy (dot-dash) processes of dusty plasma in the 
${P_d} - {V_d}$ plane for 
${\varGamma _0} = 1,{\kappa _0} = 1.5,\varGamma _0^\ast = 0.2$. In the limit of small 
${V_d}$ all the three processes become identical with 
${P_d} \propto 1/V_d^2$.
4.2. Adiabatic process
In the case where the expansion/compression time scale of dust 
${\tau _P}$ is smaller than the dust neutral collision time scale 
${\tau _{\textrm{dn}}}$, the thermal contact of the dust with neutrals is weak. In these cases the dust can regarded as isolated with constant entropy 
$(\Delta {S_d} = 0)$. This would typically be true in experiments with low neutral gas pressure (Pilch et al. Reference Pilch, Piel, Trottenberg and Koepke2007). As stated earlier, the heat exchanges involved in local plasma perturbations must be taken into account while defining thermodynamic processes in our model. Accordingly, we define the quasi-static adiabatic process involving dust as ‘the expansion/compression processes involving isolated dust that take place in an isothermal plasma background of given temperature’. Hence, in this process, while the dust is isolated and its entropy remains constant in expansion/compression, the background plasma is not and its entropy changes. The equation of state for this process can be calculated from the equation 
$0 = \Delta {U_d} + {P_d}\Delta {V_d}$. Eliminating 
${U_d}$ and 
${P_d}$ via (3.4a) we obtain the equation 
$\Delta {T_d}/{T_d} ={-} (2/3)\Delta {V_d}/{V_d}$ which can be integrated to give 
${T_d}V_d^{2/3} = \textrm{const}\textrm{.}$, (like an ideal gas). Eliminating 
${T_d}$ using this relation in the expression for 
${P_d}$ in (3.4a) finally gives the adiabatic equation of state for gaseous dust in terms one arbitrary constant C (and given 
${T_e}$ and 
${T_i}$) as
\begin{equation}{P_d} = \frac{C}{{V_d^{5/3}}} + \frac{{q_d^2N_d^2{T_e}{T_i}}}{{{q^2}n({T_e} + {T_i})V_d^2}}.\end{equation}
The adiabat/isentrope 
${P_d} = {P_d}({V_d})$ can be plotted in the 
${P_d} - {V_d}$ plane where C is determined by any point 
${P_{d0}},{V_{do}}$ (alternately 
${T_{d0}},{V_{do}}$) on the adiabat. In terms of coupling parameters 
${\varGamma _0},{\kappa _0}$, defined earlier, the equation of state can be expressed thus
\begin{equation}{\bar{P}_d} = \frac{1}{{\bar{V}_d^{5/3}}} + \frac{3}{2}\frac{{{\varGamma _0}}}{{\kappa _0^2}}\frac{1}{{\bar{V}_d^2}}.\end{equation}
In figure 1 we show a dust adiabat for 
${\varGamma _0} = 1,{\kappa _0} = 1.5,\varGamma _0^\ast = 0.2$. If the dust volume expands against external pressure from 
${\bar{V}_{d1}} \to {\bar{V}_{d2}}$, then the work done is
\begin{equation}\bar{W} = \frac{3}{2}\left( {\frac{1}{{\bar{V}_{d1}^{2/3}}} - \frac{1}{{\bar{V}_{d2}^{2/3}}}} \right) + \frac{3}{2}\frac{{{\varGamma _0}}}{{\kappa _0^2}}\left( {\frac{1}{{{{\bar{V}}_{d1}}}} - \frac{1}{{{{\bar{V}}_{d2}}}}} \right).\end{equation}
This equation shows that in expanding from 
${\bar{V}_{d1}} \to {\bar{V}_{d2}}$, dust does extra work (second term) compared with the corresponding ideal gas. The extra work is done by the ES pressure to extract requisite amount of heat from the plasma background to maintain uniform electron and ion temperatures. The changes in entropies of electrons and ions 
$\Delta {S_e},\Delta {S_i}$ can be calculated from expressions of electrons and ion entropies given in (2.3). In the simple case where 
${T_e} = {T_i} = T$, the change in the ES part of the total plasma entropy is given by 
$\Delta S = q_d^2n_d^2\Delta {V_d}/2{q^2}n$.
4.3. Constant internal energy process
In addition to the two processes described above, we can define an additional new process where the internal energy of dust 
${U_d}$ is constant. This process will take place under the condition when the dust neutral collision time scale is of the same order as the time scale of the process, i.e. 
${\tau _{\textrm{dn}}} \approx {\tau _p}$. Under this condition, the rate at which the dust exchanges heat with the dust heat bath is equal to the work done by/on the dust against the external pressure so that 
$\Delta {Q_d} = {P_d}\Delta {V_d}$ and 
$\Delta {U_d} = 0$ in (3.2). The equation of state for this process can be calculated by using (3.2), (3.4a) and the condition 
$\Delta {U_d} = 0$, and is given by
\begin{equation}{T_d} = {T_{d0}} - \frac{{2q_d^2{N_d}{T_e}{T_i}}}{{3{q^2}n{V_d}({T_e} + {T_i})}},\quad {P_d} = \frac{{{N_d}{T_{d0}}}}{{{V_d}}} + \frac{{q_d^2N_d^2{T_e}{T_i}}}{{3{q^2}V_d^2n({T_e} + {T_i})}}.\end{equation}
Thus, in the constant 
${U_d}$ process, the dust temperature increases on expansion and decreases on compression. This is an interesting and somewhat counterintuitive consequence which can be verified experimentally. However, the pressure, as usual, decreases with the volume as in the case of isothermal processes, except for a factor of 1/3 in the second term. In figure 1 we show the 
${P_d} - {V_d}$ plot for constant 
${U_d}$ process. The change in the plasma entropy, as in earlier cases, is given by 
$\Delta S = q_d^2n_d^2\Delta {V_d}/2{q^2}n$, while the change in the dust entropy is given by 
$\Delta {S_d} = (1 + 3\varGamma /2{\kappa ^2}){n_d}\Delta {V_d}$. In the large 
${\bar{V}_d}$ limit, the effects due to the mean field vanish and the constant 
${U_d}$ process becomes identical to isothermal process and 
$\Delta {S_d} = {n_d}\Delta {V_d}$.
4.4. Specific heat of dust
The specific heats of the dust at constant volume and constant pressure are defined as
\begin{equation}{C_{{V_d}}} = {\left. {\frac{{\partial {Q_d}}}{{\partial {T_d}}}} \right|_{{V_d}}},\quad {C_{{P_d}}} = {\left. {\frac{{\partial {Q_d}}}{{\partial {T_d}}}} \right|_{{P_d}}}.\end{equation}
Eliminating 
$\Delta {Q_d}$ from (2.8a,b) and using (3.2) and (3.4a,b) we get
\begin{equation}{C_{{V_d}}} = 3{N_d}/2,\quad {C_{{P_d}}} = 3{N_d}/2 + \frac{{{N_d}{T_d}}}{{{V_d}}}{\left. {\frac{{\partial {V_d}}}{{\partial {T_d}}}} \right|_{{P_d}}}.\end{equation}Eliminating the partial derivative in (4.8a,b) using (3.4a) we finally obtain the specific heat at constant pressure and the ratio of two specific heats γ as
\begin{equation}{C_{{P_d}}} = \frac{3}{2}{N_d} + \frac{{{N_d}}}{{\left( {1 + \dfrac{{3\varGamma }}{{2{\kappa^2}}}} \right)}},\quad \gamma = 1 + \frac{{2/3}}{{\left( {1 + \dfrac{{3\varGamma }}{{2{\kappa^2}}}} \right)}}.\end{equation}Equation (4.9a,b) shows that γ is no longer constant but a function of dust volume and temperature. In the ideal gas or weak coupling limit (Γ→ 0, κ → ∞), γ = 5/3.
In addition to the two types of specific heats defined in (4.7a,b), we can define a third type of specific heat, i.e. the specific heat at constant internal energy 
${C_{{U_d}}}$ defined as
\begin{equation}{C_{{U_d}}} = {\left. {\frac{{\partial {Q_d}}}{{\partial {T_d}}}} \right|_{{U_d}}} = {\left. {{P_d}\frac{{\partial {V_d}}}{{\partial {T_d}}}} \right|_{{U_d}}},\end{equation}
where we have used 
$\Delta {Q_d} = {P_d}\Delta {V_d}$ for constant 
${U_d}$. This specific heat has no ideal gas analogue, though it may have a real gas analogue. Evaluating the partial derivative and substituting for pressure from the condition in (4.6a,b), we get
\begin{equation}{C_{{U_d}}} = \frac{{3{N_d}}}{2}\left( {\frac{{2{\kappa^2}}}{{3\varGamma }} + 1} \right).\end{equation}
In the limit Γ→ 0, κ→ ∞, 
${C_{{U_d}}} \to \infty$.
4.5. Dust acoustic waves
The dust acoustic waves (DAW) are analogues of ion acoustic waves in electron–ion plasma where the inertia is due to the dust mass and the screening is due to electrons and ions. The dispersion relation of these waves can be obtained directly by using equation of states, derived above, in the fluid equations.
The total pressure of the dust 
${P_d}$, which is the sum of the kinetic and the ES pressure, drives acoustic modes in dusty plasma. The dispersion of these acoustic modes driven by the total dust pressure is governed by following set of fluid equations:
\begin{equation}{\rho _d}\frac{{\textrm{d}{\boldsymbol{v}_d}}}{{\textrm{d}t}} ={-} {\bf \nabla }({P_d}),\quad \frac{{\partial {\rho _d}}}{{\partial t}} + {\bf \nabla}\,\cdot\,{\rho _d}{\boldsymbol{v}_d} = 0,\quad {P_d} = cn_d^\gamma + \frac{{q_d^2{T_e}{T_i}}}{{{q^2}n({T_e} + {T_i})\delta }}n_d^2.\end{equation}
In these equations, 
${P_d}$ is the local dust pressure which is the sum of the dust kinetic and the ES pressure and 
${\rho _d},{v_d}$, are the local dust mass density and the fluid velocity, respectively. In (4.12c), γ = 5/3 for adiabatic processes, γ = 1 for isothermal or constant 
${U_d}$ processes, δ = 1 for adiabatic or isothermal processes and δ = 3 for constant internal energy processes. Performing standard linearization of these equations about a uniform and static equilibrium and using plane wave solutions where 
${\bf \nabla } = \textrm{i}k,\partial /\partial t ={-} \textrm{i}\omega$ in (4.12a–c), (ω and k are the frequency and the wave vector of DAW) we obtain following dispersion relation of acoustic waves:
\begin{equation}\frac{{{\omega ^2}}}{{{k^2}}} = C_{\textrm{DAW}}^2 = \left[ {\frac{{\gamma {T_d}}}{{{m_d}}} + \frac{{2q_d^2{n_d}{T_e}{T_i}}}{{{q^2}n{m_d}({T_e} + {T_i})\delta }}} \right].\end{equation}
In the low neutral gas pressure regime 
$({\nu _{\textrm{dn}}}/\omega \ll 1)$, γ = 5/3 and δ = 1 in (4.13) which corresponds to the adiabatic DAW; in the high neutral gas pressure regime 
$({\nu _{\textrm{dn}}}/\omega \gg 1)$, γ = 1 and δ = 1 which corresponds to the isothermal DAW. While in the intermediate neutral gas pressure regime 
$({\nu _{\textrm{dn}}}/\omega \approx 1)$, γ = 1, δ = 3 which corresponds to a new mode, i.e. the constant internal energy DAW.
As stated before, in DAW, the inertia is due to dust mass while the screening is due to both electrons and ions. In the usual dispersion relation of ion acoustic modes in electron–ion plasma given by 
${\omega ^2}/{k^2} = C_{\textrm{IA}}^2 = (\gamma {T_i} + {T_e})/{m_i}$, the inertia is due to ions and the screening is due to electrons. If we consider screening due to electrons in the dispersion relation of DAW given in (4.13) by taking 
${T_i} \gg {T_e}$ and the inertia due to ions by replacing dust with ions, i.e. 
${q_d} \to {q_i} = q,{m_d} \to {m_i}$, 
${n_d} \to {n_i},{T_d} \to {T_i}$, 
$2{n_i} \approx n$ (n is plasma density of electron–ion plasma) the dispersion relation of DAW in (4.13) reduces to the dispersion relation of ion acoustic given above.
In dusty plasma experiments, DAW has been observed over a wide range of frequencies and neutral gas pressures. For example, in a 3PDX device (Thomas Reference Thomas2010), DAWs in the range from 7 Hz to 120 Hz were observed. The dust particle velocity distribution function was observed to be close to Maxwellian during the wave motion (Fisher & Thomas Reference Fisher and Thomas2010; Thomas Reference Thomas2010), justifying the quasi-static assumption. The neutral pressure in this device was 72 mTorr and the corresponding dust neutral collision frequency was 
${\nu _{\textrm{dn}}} \approx 75\;\textrm{Hz}$. Then, in this experiment, waves observed in the range ω < 75 Hz will correspond to isothermal DAW, waves with ω > 75 Hz will correspond to Adiabatic DAW, while waves with ω in the range ≈75 Hz will correspond to constant internal energy DAW as discussed here. Thus it should be possible to identify the constant 
${U_d}$ DAW in this experiment which has the novel and distinctive feature that in crests dust will be cooler and trough it will be hotter.
4.6. Free expansion of dust
In dusty plasma, the hydrodynamic free expansion of Coulomb balls and dust clouds (starting from an equilibrium state) has been examined via MD simulations (Piel & Goree Reference Piel and Goree2013) and analytically using fluid equations (Ivlev Reference Ivlev2013). In experiments (Barkan & Merlino Reference Barkan and Merlino1995; Antonova et al. Reference Antonova, Du, Ivlev, Annaratonne, Hou, Kompaneets, Thomas and Morfill2012), free expansion of dust in dusty plasma has been studied in the afterglow phase of the discharge. In this phase, the confining fields are removed by switching off the anode voltage causing the background plasma to decay. The MD simulations show that at low gas pressure, the background plasma decays rapidly removing the shielding and the cloud explodes due to bare Coulomb repulsion. In contrast, at high gas pressure, the dust cloud expands gradually under the shielded Yukawa repulsion (Ivlev et al. Reference Ivlev, Kretschmer, Zuzic, Morfill, Rothermel, Thomas, Fortov, Molotkov, Nefedov and Lipaev2003; Saxena, Avinash & Sen Reference Saxena, Avinash and Sen2012) and may fission eventually (Merlino et al. Reference Merlino, Meyer, Avinash and Sen2016).
Here, in the last part of our paper, we describe the free, non-quasi-static expansion of isolated gaseous dust in a confined and steady isothermal plasma background. This is different from the free expansion of dust in an unconfined and decaying plasma background of the afterglow phase mentioned above. The non-quasi-static expansion is governed by the equation 
$\Delta {Q_d} = \Delta {U_d} + {P_{\textrm{ext}}}\Delta {V_d}$. Since the dust is isolated 
$\Delta {Q_d} = 0$ and is allowed to expand freely, 
${P_{\textrm{ext}}} = 0$. Hence the condition for free expansion is given by 
${U_d} = \textrm{const}\textrm{.}$ and the variation of temperature with volume is given by (4.6a), i.e.
\begin{equation}{T_d} = {T_{d0}} - \frac{{2q_d^2{N_d}{T_e}{T_i}}}{{3{q^2}n{V_d}({T_e} + {T_i})}}.\end{equation}
This is an irreversible process hence the entropy increases. If the dust expands freely from state 
$({T_{d1}},{V_{d1}}) \to ({T_{d2}},{V_{d2}})$ then the change in entropy is calculated by putting a reversible path between the two states. It is identical to change in entropy of the constant internal energy process and is given by 
$\Delta {S_d} = (1 + 3\varGamma /{\kappa ^2}){n_d}\Delta {V_d}$. In the ideal gas limit (Γ→ 0, κ → ∞), 
$\Delta {S_d} = {n_d}\Delta {V_d}$. Of course, as opposed to the constant internal energy process, the work done by dust in free expansion is zero. Since the dust is isolated, the free expansion must occur on time scale faster than the dust neutral time scale, i.e. 
${\tau _p} \ll {\tau _{\textrm{dn}}}$, which will be typically true at low neutral pressure. This neutral pressure condition is similar to the adiabatic process, however, the difference in the two processes is that for free expansion 
${P_{\textrm{ext}}} = 0$, while for adiabatic processes 
${P_{\textrm{ext}}} = {P_d}$. To calculate the rise in dust temperature for typical dusty plasma experimental parameters, we express (4.14) in the following dimensionless form:
\begin{equation}{\bar{T}_d} = 1 + \frac{{2{Z_d}p}}{3}\frac{{{{\bar{T}}_e}{{\bar{T}}_i}}}{{({{\bar{T}}_e} + {{\bar{T}}_i})}}(1 - {\bar{n}_d}).\end{equation}
In this equation, the dust temperature, electron and ion temperature are all normalized with initial dust temperature 
${T_{d0}}$, the dust density by the initial dust density 
${n_{d0}}$, 
${Z_d} = {q_d}/q$ and, p is given by 
$p = {Z_d}{n_{d0}}/n$. Typical parameters of dusty plasma experiments are 
${T_e} \approx 3\;\textrm{eV},{T_i} \approx 0.025\;\textrm{eV},n \approx 5 \times {10^{13}}\;{\textrm{m}^{ - 3}},{n_d} \approx 5 \times {10^9}\;{\textrm{m}^{ - 3}},{T_d} \approx 100\;\textrm{eV}$, while for micron sized dust 
${Z_d} \approx 3 \times {10^3}$ (Thomas Reference Thomas2010). Then, if we take the values of dust and plasma density given above as initial values, and if in the experiment the dust volume increases (freely) by two times so that 
${\bar{n}_d} = 1/2$, then 
${\overline T _d} = 1.15$ in (4.15), i.e. the dust temperature increases by 15 % of its initial value. In some dusty plasma experiments higher ion temperatures, 
${T_i} \approx 0.1\;\textrm{eV}$, are reported (Trottenberg et al. Reference Trottenberg, Block and Piel2006). In such cases dust temperature increases by 50 %, i.e. 
${\bar{T}_d} = 1.5$ for 
${\bar{n}_d} = 1/2$.
The reason for heating of dust particles in free expansion is the repulsive nature of the interparticle force. In free expansion of isolated dust particles, the positive potential energy decreases which is converted into thermal energy. The reverse would be true in cases of attractive force between particles. This is similar to the case of free expansion of real gases (sometimes called Joule expansion) where heating observed above inversion temperature 
${T_{\textrm{inv}}}$ is attributed to the repulsive part, and cooling below 
${T_{\textrm{inv}}}$ is attributed to the attractive part of the Lennard–Jones potential (Goussard & Roulet Reference Goussard and Roulet1993). For example, hydrogen gas with low inversion temperature (~ 200 K) shows heating on free expansion at normal temperatures. For purely repulsive potential like the present case of negatively charged dust, 
${T_{\textrm{inv}}} \to 0$, and temperature always increases in free expansion.
In the case of Joule–Thompson effect as well, where the gas is throttled through a constricted passage, the heating observed in the case 
$T \gt {T_{\textrm{inv}}}$ is due to the repulsive part of the interparticle force. In a separate paper, we will consider a Joule–Thompson-like throttling experiment for dusty plasma where dust particles are electrostatically throttled through a negatively biased mesh. This is expected to further enhance the dust temperature increase.
5. Summary and discussions
In this paper, we have proposed a thermodynamic model of dusty plasma for the case where dust is confined in a small volume within a large plasma background by external fields. The model takes into account the heat exchanged with the plasma background during the quasi-static motion of the dust. It is solved analytically in the mean field limit and various processes of gaseous phase of dust, e.g. adiabatic, isothermal and constant internal energy expansion/contraction, specific heat, dispersion of acoustic waves and free expansion of dust in the plasma background are studied.
Next, we compare our model with the HF model. The basic difference in the two models is due to the difference in the dust charge neutralization by the uniform plasma background. In the HF model where 
${V_d}/V = 1$, the dust charge is completely neutralized by the uniform plasma background. In the present model where 
${V_d}/V \ll 1$, the neutralization of dust charge by the plasma is vanishingly small. This is reflected in two different charge densities used in the HF and the present model. To show this we start with the charge density of the HF model given by (appendix B)
\begin{equation}\rho ={-} {q_d}\sum\limits_{i = 1}^{{N_d}} {\delta (\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{i}}}) + {q_d}(\overline {{n_i}} - \overline {{n_e}} ) - {\varepsilon _0}\varphi /\lambda _d^2} ,\end{equation}
where 
${\bar{n}_\alpha } = (1/V)\int {{n_\alpha }\,\textrm{d}V} = {N_\alpha }/V$ (α represents electrons and ions). Using overall charge neutrality in V, i.e. 
${q_d}{N_d} = q({N_i} - {N_e})$ we may express (5.1) in the form
\begin{equation}\rho ={-} {q_d}\sum\limits_{i = 1}^{{N_d}} {\delta (\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{i}}}) + {q_d}{n_d}\frac{{{V_d}}}{V} - {\varepsilon _0}\varphi /\lambda _d^2} .\end{equation}
From this common expression, we can obtain charge densities of both models as follows. In case 
$V = {V_d}$, we obtain the charge density of HF model. Further, using 
$\int_V {\rho \,\textrm{d}V} = 0$, we can show that the first term cancels with the second, implying thereby that the dust charge is completely neutralized (or confined) by the second term due to the uniform plasma background and the mean ES potential 
$\bar{\varphi }$ is zero; there is no need of external confinement. In contrast, in case if we take the limit of large plasma volume, i.e. 
${V_d}/V \to 0$, (for given 
${q_d}{n_d}$), then the second term in (5.2) drops out and we obtain the charge density of the present model given in (2.4). In this case the neutralization of dust charge by the uniform plasma background is asymptotically small. Further, in the dust cloud there is a non-zero mean ES potential given by 
$\bar{\varphi } ={-} {q_d}{n_d}\lambda _d^2/{\varepsilon _0}$. This negative potential tries to expel dust particles from 
${V_d}$ requiring an external field for the dust confinement. Substituting 
$\bar{\varphi } ={-} {q_d}{n_d}\lambda _d^2/{\varepsilon _0}$ in 
${P_{\textrm{ES}}} = {\varepsilon _0}{\bar{\varphi }^2}/2\lambda _d^2$ (Avinash Reference Avinash2010b) we obtain the ES part of the total dust pressure 
${P_{\textrm{ES}}}$ given by the second term in (3.4a). If we retain small terms of order 
${V_d}/V$, then 
${P_{\textrm{ES}}}$ is given by (appendix B)
\begin{equation}{P_{\textrm{ES}}} = \frac{{q_d^2N_d^2{T_e}{T_i}}}{{{q^2}nV_d^2({T_e} + {T_i})}}\left( {1 - \frac{{{V_d}}}{V}} \right).\end{equation}
From the above discussion it follows that the HF model is valid in the limit 
${V_d}/V \to 1$ while the present model is valid in the limit 
${V_d}/V \to 0$. In dusty plasma experiments 
${V_d}/V \approx {10^{ - 3}}$ to 10−5 (Barkan & Merlino Reference Barkan and Merlino1995; Trottenberg et al. Reference Trottenberg, Block and Piel2006; Pilch et al. Reference Pilch, Piel, Trottenberg and Koepke2007; Thomas Reference Thomas2010) hence neutralization of the dust charge by the plasma background is vanishingly small and an external field is used in these experiments for the dust confinement.
The expressions for the Helmholtz free energy, the dust internal energy and the dust pressure of HF model, in the mean field limit, are obtained in appendix B and are given by
\begin{gather}{F_d} = {T_d}{N_d}[\ln ({n_d}\varLambda _d^3) - 1] + \sum\limits_\alpha {{T_\alpha }{N_\alpha }[(\ln {n_0}\varLambda _\alpha ^3) - 1]} ,\end{gather}
\begin{gather}{U_d} = \tfrac{3}{2}({N_d}{T_d}) + \sum\limits_\alpha {{T_\alpha }{N_\alpha }[(\ln {n_0}\varLambda _\alpha ^3) - 1]} ,\end{gather}
\begin{gather}{P_d} = \frac{{{N_d}{T_d}}}{{{V_d}}} = {n_d}{T_d}.\end{gather}Clearly, in the HF model, the ES contributions are zero in the mean field limit, i.e. the mean field limit contains only thermal terms; ES contributions arise solely due to finite dust correlation effects. In contrast, in our model given in (2.3) and (2.4), thermodynamic parameters have a dominant ES contribution even in absence of dust correlations (or in the mean field limit).
The ES part of the dust pressure 
${P_{\textrm{ES}}}$ has been measured in experiments (Fisher et al. Reference Fisher, Avinash, Thomas, Merlino and Gupta2013; Williams Reference Williams2019) where, especially for smaller clouds, it is found to be substantially greater than the thermal dust pressure predicted by the HF model in (5.6). In fact, it is of the same order as the pressure predicted by our model in (3.4a) which has substantial ES contribution (Fisher et al. Reference Fisher, Avinash, Thomas, Merlino and Gupta2013). Additionally, in the HF model the dust density is assumed to be proportional to average plasma density for which there is no experimental justification. In experiments, dust and plasma density can be varied independently (Barkan & Merlino Reference Barkan and Merlino1995; Thomas Reference Thomas2010) as in our model. In table 1 we show the comparison of the two models.
Table 1. Comparison of the present model with the HF model.
Finite dust correlation effects are contained in the double summation term of the effective dust internal energy expression given in (2.9). In the present paper we have appropriately approximated the double summation term to construct the thermodynamic limit of dusty plasma. A theory of finite correlation effects, e.g. melting/freezing phase transitions, sound propagation in correlated medium etc. can be constructed analytically by suitable approximation of the double summation term or by calculating this term via MD simulations. This will be the subject of future publications.
Acknowledgements
The author is grateful to one of the referees for thoughtful comments. This work is dedicated to the memory of P.K. Kaw, my teacher and mentor, who greatly appreciated and advocated the use of simple thermodynamic arguments in plasma physics pioneered by early workers in the field.
Editor Edward Thomas thanks the referees for their advice in evaluating this article.
Declaration of interests
The author reports no conflict of interest.
Appendix A. Derivation of internal energy
The expression for the composite internal energy of the system of plasma and dust is given by
\begin{equation}U = \frac{3}{2}({N_e}{T_e} + {N_i}{T_i} + {N_d}{T_d}) + \frac{1}{2}\int {\rho \varphi \,\textrm{d}r} - \frac{{q_d^2}}{{8{\rm \pi} {\varepsilon _0}}}\sum\limits_{j = 1}^{{N_d}} {\int {\frac{{\delta (\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{j}}})}}{{|\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{j}}}|}}} \,\textrm{d}r,}\end{equation}where ρ and φ are given by
\begin{equation}\rho ={-} {q_s}\sum\limits_{i = 1}^{{N_d}} {\delta (\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{i}}}) - {\varepsilon _0}\varphi /\lambda _d^2} ,\quad \varphi ={-} \frac{{{q_d}}}{{8{\rm \pi} {\varepsilon _0}}}\sum\limits_j {\frac{{{\rm exp} ( - |\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{j}}}|/{\lambda _d})}}{{|\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{j}}}|}}}. \end{equation}Eliminating ρ and φ through (A 2) in the second term of (A 1), and performing the integration with delta function (after subtracting the singular term) we obtain
\begin{align} U &= \dfrac{3}{2}({N_d}{T_d} + {N_i}{T_i} + {N_e}{T_e}) - \left( {\dfrac{{q_d^2{N_d}{\kappa_d}}}{{8{\rm \pi} {\varepsilon_0}}}} \right) + \left[ {\dfrac{{q_d^2}}{{8{\rm \pi} {\varepsilon_9}}}\sum\limits_i {\sum\limits_{j \ne i} {\left( {\dfrac{{{\rm exp} ( - {\kappa_d}|{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|)}}{{|{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|}}} \right)} } } \right]\nonumber\\ & \quad - \dfrac{{{\varepsilon _0}\kappa _d^2}}{2}\int {{\varphi ^2}{\textrm{d}^3}r} . \end{align}To evaluate the last integral in (A 3), we eliminate φ in the φ 2 integral in (A 3) by (A 2b), which gives
\begin{equation}\frac{{{\varepsilon _0}\kappa _d^2}}{2}\int {{\varphi ^2}{\textrm{d}^3}r} = \frac{{{\varepsilon _0}\kappa _d^2}}{2}{\left( {\frac{{{q_d}}}{{4{\rm \pi} {\varepsilon_0}}}} \right)^2}\int {\sum\limits_i {\frac{{{\rm exp} ( - {\kappa _d}|\boldsymbol{r} - {\boldsymbol{r}_{{\boldsymbol{i}^{}}}}|)}}{{|\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{i}}}|}}} \sum\limits_j {\frac{{{\rm exp} ( - {\kappa _d}|\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{j}}}|)}}{{|\boldsymbol{r} - {\boldsymbol{r}_{{\boldsymbol{j}^{}}}}|}}\,{\textrm{d}^3}r.} }\end{equation}The integral is evaluated using spherical polar coordinates (Avinash Reference Avinash2010a) to give
\begin{equation}\sum\limits_i {\sum\limits_j {\int {\frac{{\rm exp} ( - {\kappa _d}|\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{i}}}|)}{|\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{i}}}|}\frac{{\rm exp} ( - \kappa _d|\boldsymbol{r} - \boldsymbol{r}_{\boldsymbol{j}}|)}{|\boldsymbol{r} - \boldsymbol{r}_{\boldsymbol{j}}|}\,{\textrm{d}^3}r} } } = \frac{2{\rm \pi}}{\kappa _d}\sum\limits_i {\sum\limits_j{^{{'}}{\rm exp} ( - {\kappa _d}|{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|)} } .\end{equation}Eliminating the integral in (A 3) using (A 4) and (A 5) gives the (2.6) of the main text.
To derive (2.7) of the text, we start with expression of electron and ion entropy in (2.3) of the text and eliminate the electron and ion densities through linearized Boltzmann response. The resulting expression is expanded in powers of 
$q\varphi /{T_e} \approx ,q\varphi /{T_i} \lt 1$ and terms proportional to φ 2 are retained to obtain following expression (Hamaguchi & Farouki Reference Hamaguchi and Farouki1994):
\begin{equation}{T_e}{S_e} + {T_i}{S_i} = \frac{3}{2}\sum\limits_\alpha {{N_\alpha }{T_\alpha }} - \sum\limits_\alpha {{T_\alpha }{N_\alpha }(\ln n\varLambda _\alpha ^3 - 1) - \frac{{{\varepsilon _0}\kappa _d^2}}{2}\int {{\varphi ^2}\,{\textrm{d}^3}r} } .\end{equation}Eliminating the last integral in (A 6) through (A 4) and (A 5) we obtain (2.7) of the main text.
Appendix B. Homogenous limit of the Hamaguchi–Faouki model
In the HF model the background electron–ion plasma and dust particles occupy the same volume V.
The linearized Boltzmann response of HF model is given by 
${n_\alpha } = {\bar{n}_\alpha }(1 - {q_\alpha }\varphi /{T_\alpha })$ where 
${\bar{n}_\alpha } = (1/V)\int {{n_\alpha }\,\textrm{d}V} = ({N_\alpha }/V),\int {\varphi \,\textrm{d}V} = 0$ and α denotes ions or electrons. With this linearized response, and using the overall charge neutrality in V given by 
${q_d}{n_d} = q({\bar{n}_i} - {\bar{n}_e})$, the net charge density is given by
\begin{equation}\rho ={-} {q_d}\sum\limits_{i = 1}^{{N_d}} {\delta (\boldsymbol{r} - {\boldsymbol{r}_{\boldsymbol{i}}}) + {q_d}{n_d} - {\varepsilon _0}\varphi /\lambda _d^2} .\end{equation}
The corresponding expression for the dust internal energy 
${U_d}$ of the HF model, without periodic boundary condition terms, is given by eliminating ρ from (B 1) in (2.2) of the text, which gives
\begin{align}{U_d} & = \dfrac{3}{2}({N_d}{T_d}) + \dfrac{{q_d^2}}{{8{\rm \pi} {\varepsilon _0}}}\sum\limits_i {\sum\limits_j {\dfrac{{{\rm exp} ( - |{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|/{\lambda _d})}}{{|{\boldsymbol{r}_{\boldsymbol{i}}} - {\boldsymbol{r}_{\boldsymbol{j}}}|}}} } + \sum\limits_\alpha {{T_\alpha }{N_\alpha }[(\ln {n_0}\varLambda _\alpha ^3) - 1]} \nonumber\\ &\quad - \dfrac{{q_d^2N_d^2{T_e}{T_i}}}{{{q^2}{V_d}({T_e} + {T_i})}} - \dfrac{{{N_d}q_d^2{\kappa _d}}}{{8{\rm \pi} {\varepsilon _0}}}. \end{align}
The second last term (which is absent in our model) arises due to 
${q_d}{n_d}$ in (B 1). It corresponds to a cohesive field in the uniform plasma background which, as stated earlier, neutralizes or confines the negative dust charge. If now we take the homogenous limit in (B 2), then the second term cancels with the second last term. The corresponding expressions for 
${U_d}$ and the Helmholtz free energy in the homogenous limit, which now contain only thermal terms, are
\begin{gather}{U_d} = \tfrac{3}{2}({N_d}{T_d}) + \sum\limits_\alpha {{T_\alpha }{N_\alpha }[(\ln {n_0}\varLambda _\alpha ^3) - 1]} ,\end{gather}
\begin{gather}{F_d} = {T_d}{N_d}[\ln ({n_d}\varLambda _d^3) - 1] + \sum\limits_\alpha {{T_\alpha }{N_\alpha }[(\ln {n_0}\varLambda _\alpha ^3) - 1]} ,\end{gather}where α denotes electrons and ions. The dust pressure is given by
\begin{equation}{P_d} = \frac{{{N_d}{T_d}}}{{{V_d}}}.\end{equation}
To obtain 
${P_{\textrm{ES}}}$ with small corrections of order 
${V_d}/V \ll 1$, we calculate 
$\bar{\varphi }$ from (5.2), retaining the 
${V_d}/V$ corrections (middle term) to give
\begin{equation}\bar{\varphi } ={-} \frac{{{q_d}{n_d}\lambda _d^2}}{{{\varepsilon _0}}}\left( {1 - \frac{{{V_d}}}{V}} \right).\end{equation}
It has been shown earlier that 
${P_{\textrm{ES}}}$ can be obtained directly from the expression 
${P_{\textrm{ES}}} = \int {{q_d}{n_d}\,\textrm{d}\bar{\varphi }}$ (Avinash Reference Avinash2010b). Eliminating 
${q_d}{n_d}$ in the integral from (B 6) we obtain
\begin{equation}{P_{\textrm{ES}}} = \frac{{q_d^2N_d^2{T_e}{T_i}}}{{{q^2}nV_d^2({T_e} + {T_i})}}\left( {1 - \frac{{{V_d}}}{V}} \right),\end{equation}
where we have eliminated Debye length using 
$1/\lambda _d^2 = ({q^2}{n_0}/{\varepsilon _0})(1/{T_e} + 1/{T_i})$.
