You are an expert at summarizing long articles. Proceed to summarize the following text: 
gas and dust mixtures are ubiquitously present in nature , so a correct numerical prescription of its evolution turns out to be of the uttermost importance . in many astrophysical applications , dust can be described as a set of particles immersed in a fluid phase ( gas ) . mathematically , such a system can be described using the @xcite notation , by the following set of equations @xmath0 @xmath1 @xmath2 @xmath3 @xmath4 where @xmath5 and @xmath6 are the dust particles number density and mass respectively , @xmath7 is the gas density , @xmath8 and @xmath9 are the dust and gas velocities , @xmath10 is the gas thermal energy , @xmath11 is the drag coefficient for a single particle , @xmath12 represents the gas pressure , and @xmath13 stands for any external forces , like gravity or radiation pressure . note than in equation [ eu2 ] , the external force per unit volume @xmath14 is required for the gas . @xmath15 is the lagrangian derivative , and its specific form will be discussed in section 2 . the effects of forces related to the intrinsic volume of the dust particles have been ignored , since in normal astrophysical applications they become negligible .    in the present work , we will concentrate on the study of drag forces . the form of the drag force of gas on a single dust grain may vary considerably as a function of the grain and gas properties @xcite . if the mean free path of the gas molecules is bigger than the dust particle radius s ( assuming spherical grains ) , the expression of the drag coefficient on a single dust grain becomes        @xmath16    where @xmath17 is the velocity of the gas molecules due to thermal motion , @xmath18 the gas temperature , @xmath19 is the mean molecular weight and @xmath20 is the atomic mass of hydrogen . if on the contrary , the mean free path of the gas molecules is smaller than the dust particle radius , the expression of the drag force on a single dust particle becomes @xmath21 where the dimensionless coefficient @xmath22 will be given by @xcite @xmath23 @xmath24 @xmath25 where @xmath26 is the reynolds number and @xmath27 is the molecular viscosity of the gas . under certain circumstances ( typically for small dust grain sizes ) , the acceleration experienced by the dust can become very large , leading to very short stopping times . the occurrence of such short stopping times may become , under certain circumstances , a very severe problem in the numerical simulation of dust and gas mixtures . in protoplanetary disks , for example , the typical range of body sizes spreads from micron - sized dust grains , up to kilometre - sized planetesimals . consequently , the ranges of dust - gas coupling intensities and stopping times will be large , leading to a large range of dynamical time scales .    the first attempt to study gas and dust mixtures in the framework of the sph method was developed by @xcite , and was subsequently improved by @xcite by the inclusion of an implicit time - integration scheme . the main problem with the method was its incapacity to guarantee a convergent solution under certain circumstances . @xcite proposed a variation of @xcite method . despite being capable of providing stable and convergent solutions , their method still suffers three main difficulties , intrinsic to any two fluid approach : ( i ) an inclination to produce artificial dust clumps whenever the dust is concentrated below the gas resolution , due to the pressureless nature of the dust component , ( ii ) the necessity of a very high spatial resolution , in the high drag regime in order to avoid overdissipation , and ( iii ) the necessity of a very high number of iterations in the implicit time integration scheme , or a very high number of time - steps in the explicit scheme , for the high drag regime . more recently , a new one fluid approach has been proposed by the same authors @xcite . in this new approach , both fluids are evolved as a single fluid by using the barycentric velocity as a common reference frame . through this approach , most of the aforementioned problems are avoided . however , in its present state , the one fluid method struggles with the low drag regime when dust and gas are not well described as a mixture and the velocity field should be multi - valued @xcite , whereas a two - fluid method handles this situation with ease . in this paper , a new two - fluid sph method will be investigated in order to solve the third of the aforementioned problems . a simple semi - analytical model is proposed , in order to approximate the time evolution of the dust component , and thus avoid the need for a numerical integration of its time evolution . special attention will also be paid to the impact of overdissipation in the method . in particular , it will be shown that the method is much better at resolving dusty shocks in the limit of short stopping times than other explicit or implicit two - fluid sph methods . whenever possible , an estimation of the resolution requirements of the method will be provided . this paper is organized as follows . in section 2 , the possibility of imposing an analytical decay model as an approximate solution for the small dusty grains evolution will be discussed . in section 3 , a series of numerical tests will be presented in order to compare the accuracy of the present method with more traditional approaches . finally , in section 4 , we will draw our conclusions .      + as mentioned in the introduction , the objective of the present work is to avoid the need for a full numerical integration of the velocity evolution of small dust grains , whenever the stopping time becomes prohibitively short . in order to do so , one could try to estimate the total change in velocity of a dust particle , after having interacted through drag with the gas , for a certain time @xmath28 . as seen in the introduction , if we concentrate exclusively on the drag interaction , the equations of motion for the time - evolution of an arbitrary pair of dust and gas fluid elements ( represented in a two - fluid sph method by a pair of particles located at positions @xmath29 and @xmath30 ) are @xmath31 @xmath32 @xmath33 where @xmath34 is the volume density of the dust component , @xmath35 and we consider the epstein regime where we have defined @xmath36 . in the present work , the adopted evolutionary equations for the dust and gas components are @xmath37 @xmath38 @xmath39 where @xmath40 equations [ vd_sol ] , [ vg_sol ] and [ ug_sol ] will constitute an approximate solution for the equations of motion , as long as dust and gas densities can be considered as approximately constant along the integration time - step @xmath28 , since @xmath41 @xmath42 @xmath43 the main attractive of equations [ vd_sol ] , [ vg_sol ] and [ ug_sol ] is that they can be used to approximately describe both strong and weak drag regimes . if @xmath44 , equations [ vd_sol ] to [ ug_sol ] become @xmath45 @xmath46 @xmath47 whereas if @xmath48 , equations [ vd_sol ] , [ vg_sol ] and [ ug_sol ] simply read @xmath49 @xmath50 @xmath51 which is the expected solution for the equations of motion of a strongly coupled dust and gas mixture . another attractive feature of equations [ vd_sol ] to [ ug_sol ] is that they naturally incorporate , due to its fully lagrangian nature , perfect advection into the numerical scheme . if one calculates the time evolution of the relative velocity between dust and gas , in the dust frame , one gets @xmath52 so , as long as the velocity evolution for each phase is calculated by using the local acceleration on each frame , the extra term related to the differential velocity of the frames , will be naturally included into the scheme . this property , although not completely intuitive , can clearly be seen if one considers the case of ballistic particles and a gas that do not interact at all ( something that sph can treat very easily ) . the lagrangian equations that describe the evolution of such a system , which are the equations that would be solved by an sph implementation , are @xmath53 @xmath54 if one now calculates the time variation of the relative velocity between the phases as in equation [ lagdg ] , we obtain @xmath55 where the second term on the right hand side just reflects that we have had to choose between the dust and the gas when defining our lagrangian derivative . it is not an extra term that needs to be implemented . note that in some recent sph one - fluid prescriptions , the extra advection terms do need to be explicitly calculated ( e.g. equation 14 of @xcite ) . this characteristic should be clearly considered as an advantage of our method . the key to our two fluid method for modelling a dusty gas is that we now operator split the differential equations that describe the evolution of gas and dust , so that we solve everything except the drag term using standard explicit integration methods , and subsequently modify the resulting velocities by applying the drag term separately . for example , to include gas pressure and drag forces between the dust and the gas , we first use standard explicit sph integration to apply @xmath56 @xmath57 @xmath58 and then , we apply equations [ vd_sol ] to [ ug_sol ] to the obtained intermediate velocities and thermal energy @xmath59 , @xmath60 and @xmath61 @xmath62 @xmath63 @xmath64      +    in order to apply equations [ vd_sol2 ] , [ vg_sol2 ] , and [ ug_sol2 ] in the sph two fluid approach , the gas and dust elements are discretized into a set of mass elements , often called particles . any continuous quantity will be thus reconstructed by means of an interpolation method @xmath65 @xmath66 where @xmath67 is the mass of each sph particle , @xmath68 is the smoothing length of each sph particle , and @xmath69 is the interpolating function , called the kernel ( see for example * ? ? ? * ) . in general , in the two fluid scheme , the value of the gas velocity at a dust location ( and vice versa ) will be unknown , so in equations [ vd_sol2 ] , [ vg_sol2 ] , and [ ug_sol2 ] the use of equations [ interp1 ] and [ interp2 ] will be necessary . in particular , using the @xmath70 index to refer to dust particles , @xmath71 to gas particles , and @xmath72 to the neighbours of opposite type , we can evaluate the difference between the dust and gas velocities as @xmath73 @xmath74 where @xmath75 , @xmath76 , @xmath77 , and @xmath78 . by using sph interpolation , equations [ vd_sol2 ] , [ vg_sol2 ] and [ ug_sol2 ] can be discretized @xmath79 @xmath80 @xmath81 where @xmath82 and @xmath83 are normalisation factors ( see below ) . physically , sph particles should be understood as finite mass elements of each one of the components . in particular , sph dust particles should be interpreted as homogeneous ensembles of dust particles of radius @xmath84 , intrinsic mass @xmath6 , and number density @xmath5 . therefore , for each sph dust particle one can assign a volume density @xmath85 which will represent the total dust mass contained within the volume of the sph dust particle ( determined by its kernel support radius ) . smoothing lengths and volume densities for both components can be calculated by a standard iterative sph manner , solving @xmath86 through a newton - raphson method @xcite , where @xmath87 for the standard cubic spline kernel , and the dust and gas densities are given by @xmath88 @xmath89 this procedure is equivalent to solving the continuity equations [ cont1 ] and [ cont2 ] . sph particle masses will be assigned by dividing the total mass of each component present in the simulation , by the number of particles of the component .    in order to calculate the dust - to - gas ratio at a given dust particle location , we estimate the gas and dust mass fraction contained within the interpolation sphere of the sph dust particle . that is , we take @xmath90 this prescription is chosen due to its greater stability , in comparison with the simpler dust and gas densities quotient . we have found that , whenever discontinuities are present in the computational domain ( for example in the shock tube test ) , the fluctuations in the dust density can lead to high stopping time fluctuations if the ratio @xmath91 is used directly in equation [ ts ] . if equation [ dust - to - gas ] is used , because the mass of the sph dust particle is constant , the fluctuations are avoided . furthermore , this approach allows us to calculate dust evolution even with a very low number of sph dust particles , since it does not rely on the validity of the fluid approximation for the dust component . also , and in order to minimize fluctuations if a low number of neighbours is present , a normalization factor @xmath82 has also been included in the sph dust summation @xcite , equal to @xmath92 due to the symmetric structure of equations [ vd_sph1 ] and [ vg_sph1 ] linear momentum is preserved during the interaction , and a projection of the relative velocity along the line of sight of the particlesis introduced in order to guarantee angular momentum conservation @xcite . a normalization factor @xmath27 , equal to the number of the spatial dimensions of the system , is necessary to guarantee the equivalence of the projection method with equations [ vd_sol ] , [ vg_sol ] and [ ug_sol ] up to a second order approximation ( see @xcite for an excellent discussion ) . for the same reason , energy can also be shown to be conserved . the kinetic energy of the mixture , at @xmath93 , will be expressible as @xmath94 so the change in kinetic energy will be @xmath95 then , by assuming that @xmath96 one finds , after introducing the sph summations @xmath97 the method has been tested with two different integrators , a second order runge - kutta fehlberg @xcite , and a second order predictor - corrector @xcite . the obtained results with the two integrators have been equivalent in all cases , except in the sound wave test ( section 3.2 ) where the runge - kutta scheme leads to a poorer energy and momentum conservation . please , see appendices a and b for a detailed explanation of both integration methods .      to investigate the stability of the numerical scheme , equations [ vd_sol ] and [ vg_sol ] may be written in the following form @xmath98 @xmath99 where @xmath100 as can be seen , equations [ vd_stab ] and [ vg_stab ] can be interpreted as a forward euler method , where velocity is evolved with respect to @xmath101 instead of time . following @xcite , a von newmann analysis can be done . if the dust and gas components are perturbed with a monochromatic plane wave @xmath102 [ vd_vn ] @xmath103 equations [ vd_stab ] and [ vg_stab ] may be written as the following linear system @xmath104 the corresponding two eigenvalues of the system are @xmath105 and the system will remain numerically stable ( @xmath106 ) whenever @xmath107 which will always occur , given the definition of @xmath101 , except in the limit @xmath108 . in this case , equation [ cour1 ] will act as a courant - like condition . in order to keep stability , it will be enough to decrease the integration @xmath101-step by a factor of 2 , and evolve the system of equations [ vd_stab ] and [ vg_stab ] in two steps . because of the existing linear relation between the velocity and @xmath101 , the accuracy of the solution will not be affected by the number of steps performed , like in an ordinary explicit integration scheme . the present method possesses two very different regimes , depending on the ratio between the gas integration time - scale @xmath28 , and the dust stopping time @xmath109 . if @xmath110 ( when an explicit integration could be used ) @xmath111 and equations [ vd_sph1 ] , [ vg_sph1 ] and [ ug_sph1 ] thus become @xmath112 @xmath113 @xmath114 in this limit , the stability condition [ cour1 ] becomes @xmath115 coinciding with the courant condition of an explicit integration as shown by @xcite . therefore , equations [ vd_sph3 ] to [ ug_sph3 ] will be equivalent to an explicit sph two fluid method @xcite as long as very sharp density gradients are absent from the gas component . to quantify the errors produced by this approximation , the behaviour of the algorithm in the presence of strong density gradients ( shocks ) will be tested in section 3.3 . it is also interesting to check whether the method can reproduce the properties of the dust and gas mixture in the so - called _ terminal velocity approximation _ ( see for example @xcite and references therein ) . when dust and gas are strongly coupled , the dust reaches a constant relative velocity with respect to the gas , which is small but still finite . such a relative velocity is proportional to the pressure gradient and the stopping time @xmath109 . one can see this by using @xmath121w(r_{\rm kj } ) \end{aligned}\ ] ] where we have introduced @xmath122 to simplify the notation . now , in order to reach the terminal velocity , pressure gradient and drag forces must balance each other , leading to @xmath123w(r_{\rm kj})\\ & = \sum_{\rm k}^{dust}\frac{m_{\rm k}}{\rho_{\rm k}}\left[-\frac{\nabla p_{\rm j}}{\rho_{\rm j}}+\frac{1}{t_{\rm s}}(\textbf{v}_{\rm kj}\cdot\hat{\textbf{r}}_{\rm kj})\hat{\textbf{r}}_{\rm kj } \right]w(r_{\rm kj})\end{aligned}\ ] ] which is simply the sph equivalent of    @xmath124      the procedure followed in section 2.1 can be extended to the non - linear drag regimes as long as an approximate analytic solution can be found for the time evolution of the dust grains . for example , in a full non - linear regime ( equations [ kst ] and [ cd3 ] ) , a procedure analogous to the one in section 2.1 can be followed . in such a regime , the equations of motion of the dust and gas components can always be expressed as    @xmath125    @xmath126    where @xmath127 . in this case , the chosen equations for a pair of arbitrary dust and gas fluid elements located at points @xmath29 and @xmath30 are @xmath128 @xmath129 @xmath130 where @xmath131         by using sph discretization , equations [ vd_sol3 ] , [ vg_sol3 ] and [ ug_sol3 ] become @xmath132 @xmath133 @xmath134 where in this case , an additional sph summation is necessary to calculate @xmath135 , since it depends on the relative velocity of the components at the dust particle location . to perform most of the numerical tests , the dragging algorithm was implemented in a purpose - built sph code . the code included self - consistent @xmath136 and @xmath137 calculation , grad - h terms @xcite , and riemann solver - like artificial viscosity with thermal conductivity whenever neeeded @xcite . to perform the sedov test , the dragging algorithm was implemented into a well tested three - dimensional sph code . for the sake of conciseness , the exact details of the sph code will not be presented here , but the interested reader is referred to @xcite .          the dustybox test @xcite was performed in order to prove the capacity of the method to reproduce the expected asymptotic behaviour of the drag force . a set of @xmath138 dust and @xmath138 gas particles with homogeneous densities @xmath7 and @xmath85 are placed in a periodic box with an initial velocity @xmath139 and @xmath140 . in order to construct the initial model , particles are evenly distributed along a cubic lattice with @xmath141 . the dust lattice is shifted , with respect to the gas , by half of the gas particles separation in each direction . the mass of each sph particle is equal to @xmath142 where @xmath143 is the computational domain volume , and @xmath144 the number of particles in each phase . an isothermal equation of state is adopted ( @xmath145 ) , and in this case no artificial viscosity is used . the physical units of the problem are chosen such that @xmath146  g  @xmath147 , @xmath148  g  @xmath147 , @xmath149  cm  s@xmath150 . these are the appropriate conditions for a dust particle at the mid - plane of a protoplanetary disk at 1 au from the central star ( see for example @xcite ) . the computational domain comprises a total volume of 1 cubic au , and the total mass of gas inside the domain is @xmath151  g. the integration time - step @xmath28 is calculated by finding the minimum value , for all gas particles , of @xmath152 and @xmath153 where @xmath137 is the sph particle smoothing length and @xmath154 is the gas particle acceleration . since the pressure gradient is zero , the exact solution for equation ( 4 ) is easy to find in this case @xmath155 allowing a direct comparison of the results obtained . in fig . [ fig : test1a ] the time evolution of the velocity of a single sph dust particle is presented for several different values of the dust grain size @xmath84 . left figure corresponds to a case with @xmath156 while right figure corresponds to a case with @xmath157 . as can be seen , irrespectively of the dust grain size , the correct terminal velocity between gas and dust components is reached in all cases . whenever the gas integration time step becomes smaller than the dust stopping time , the algorithm is capable of following the velocity decay of the dust component towards its limiting velocity . if the dust stopping time becomes much smaller than the gas integration time step , the algorithm simply tries to put both components on their barycentric velocity , right from the start . if under any circumstance , resolving the velocity decay becomes essential , one can always artificially decrease the gas integration time - step by reducing the gas courant time condition ( equation [ courant0a ] ) by an arbitrary factor . in fig . [ fig : test1b ] time evolution of the dust component velocity in the three - dimensional @xmath157 , @xmath158  mm case is shown , for different values of the gas integration time step . as can be seen , since the stopping time is much shorter than the gas integration time - step ( @xmath159  yrs ) , an artificially reduced gas integration time - step is needed to start resolving the dust component velocity decay .    in order to better appreciate the precision of the adopted approximation , in fig . [ fig : test1c ] , the relative errors in the @xmath157 test with @xmath158  m ( left plot ) and @xmath158  mm ( right plot ) are presented , for several different kernels . we test the standard @xmath160 cubic spline kernel @xcite , the @xmath161 quintic spline kernel , and the double hump cubic kernel @xcite . as can be seen , the correct terminal velocity is obtained , irrespectively of the used kernel , with very high precision ( the relative error between the numerical and analytical results is @xmath162 ) . the greatest departures from the analytical solution are obtained during the velocity decay phase . in this phase , only the double hump kernel keeps errors under acceptable limits . a similar result was also found by @xcite in their study . one can also see from the right plot of fig . [ fig : test1c ] , that the normalization condition ( equation [ norm ] ) , helps to reduce the errors further . if the double hump kernel is used in conjunction with the normalization condition , the maximum relative error is @xmath163 in all the tested cases . this is a very important results , since it shows that the trajectories of dust particles of arbitrary size can be accurately calculated , in a protoplanetary - like environment , without the need for excessive resolution . the second test performed was the study of the propagation of a sound wave in a dust - gas mixture in a constant drag regime , also known as the dustywave test @xcite . this can be done by setting the drag coefficients on a single grain , in the epstein regime , to be equal to @xmath164 as a consequence , the equations of motion for the dust and gas components become @xmath165 @xmath166 where we have introduced the drag coefficient per unit volume @xmath167 . this is a particularly interesting problem , because as in the previous case , it possesses an analytic solution @xcite . to set up the test , an ensemble of dust and gas particles with homogeneous densities @xmath85 and @xmath7 are evenly distributed over a periodic one - dimensional domain @xmath168 particle masses are assigned in the same way as in the previous section . no artificial viscosity is used in this case in order to avoid introducing non - physical energy dissipation in the test . the integration time - step @xmath28 is again calculated by finding the minimum value given by equations [ courant0a ] and [ courant0b ] . an isothermal equation of state @xmath169 with @xmath170 is used in this case . in order to create the waves , a sinusoidal perturbation is introduced for each particle , both in position and velocity @xmath171 @xmath172 where @xmath173 is the original position of each particle , and @xmath174 , so that the velocity perturbation of the wave is @xmath175 . the spatial perturbation @xmath176 will be different for every resolution , and is selected in each case so that the density perturbation of the wave is always @xmath177 . after introducing the perturbation , the propagation of the resulting sound wave within the domain is followed . as previously mentioned , we used the predictor - corrector integrator for this test , since it gives better long - term energy and momentum conservation .      +    in fig . [ fig : test2a ] , four different snapshots of the time evolution of the sound wave velocity in a @xmath178 , @xmath179 case with @xmath180 ( @xmath181 ) ( left panel ) , and @xmath182 ( @xmath183 ) ( right panel ) , are presented . in fig . [ fig : test2b ] , four different snapshots of the time evolution of the sound wave velocity in a @xmath184 , @xmath185 case with @xmath186 ( @xmath187 ) ( left panel ) , and @xmath188 ( @xmath183 ) ( right panel ) , are presented . as can be seen , a good agreement with the analytical solutions has been obtained in both cases . in order to quantify the deviation from the analytical solution several error norms are calculated ( see figure captions ) @xmath189 @xmath190^{1/2},\ ] ] @xmath191 where @xmath192 is the maximum value of the exact solution in the plotted region , @xmath193 is the analytical solution for the @xmath70-th point of the plot , and @xmath144 is the number of plotted points ( http://users.monash.edu.au/dprice/@xmath194splash/userguide/ ) . as previously mentioned , one of the most important characteristics of dust and gas mixtures is that the local sound speed modification as a function of the dust / gas fraction ( equation [ cs ] ) . since the analytical solutions seen in fig . [ fig : test2a ] and  [ fig : test2b ] take into account such a modification , the test confirms the capacity of the algorithm to reproduce this feature of dust / gas mixtures . the results of figs . [ fig : test2a ] and  [ fig : test2b ] also confirm that , whenever the amount of momentum transferred between the phases is small compared with the total momentum of the gas , an arbitrarily low number of dust particles can be used . both in the lower drag case of fig . [ fig : test2a ] , and in fig . [ fig : test2b ] only 8 dust particles per wavelength are necessary to obtain reasonable results . unfortunately , as can be seen in fig . [ fig : test2c ] , a certain excess of energy dissipation by drag in the high drag regime becomes unavoidable . this is not a new phenomena and was already found by @xcite in their simulations . the sph two fluid scheme needs a minimum resolution ( @xmath195 ) in order to correctly resolve the small position and velocity differences between the dust and gas phases , otherwise overdissipation becomes unavoidable . however , because our method treats only the gas as a fluid , and not the dust , this resolution criterion must only be satisfied by the gas component , not the dust . in fig . [ fig : test2c ] , the effect of particle resolution is investigated in the high drag regime for two different dust - to - gas ratios . from top to bottom a total of 256 , 128 , 64 and 32 gas particles have been used . in the left panel of fig . [ fig : test2c ] a @xmath156 case with @xmath183 is presented . in this case , equal numbers of gas and dust particles have been used . in the right panel of fig . [ fig : test2c ] a @xmath157 case with @xmath183 is presented . in this case only 8 dust particles have been used . as can be seen , and in complete agreement with the minimum resolution condition , only the ones with a minimum number of @xmath196 gas and dust particles is capable of matching the expected solution in the @xmath156 case . however , in the @xmath157 case , overdissipation effects become much less dramatic , even with a very low gas and dust particle resolution . this is important , since most astrophysical applications have low dust - to - gas ratios . it is also important to note that the present method is less dissipative than the previous ones , because we need to perform many fewer integration time - steps , in order to evolve the simulation to a given time . in the present method , the interpolation error is only committed once per gas integration time step , in contrast with explicit or implicit methods where the error can be committed hundreds or thousands of times per gas integration time step . in fact , in earlier versions of the present method an iterative procedure was tried in order to achieve a higher precision in the final relative velocities between the components , but instead it resulted in a degree of overdissipation comparable to the one using a standard integration method . the next two tests are the shock tube test , and the sedov blast test @xcite . they were both performed in order to test the behaviour of the scheme , in the presence of strong density and pressure gradients . since equation ( 16 ) will only be valid as long as no big changes in the density or the pressure gradient occur during the integration time - step , these experiments are critical to prove the usefulness of the method . in these experiments , thermal energy plays an essential role in the evolution of the system , so this time an adiabatic equation of state with @xmath197 and @xmath198 is used . also , to correctly model the shocks , @xcite artificial viscosity is used in both cases with coefficients @xmath199 and @xmath200 for the thermal conduction parameter . the signal velocities are @xmath201 and latexmath:[$\textbf{v}_{\rm sig , u}=    time - step @xmath28 is calculated in both cases by finding the minimum value for all gas particles between @xmath203 and @xmath204      +      where @xmath154 is the sph gas particle acceleration . note that if additional forces affecting both phases ( like radiation pressure for example ) were introduced in the simulation , condition [ dt_a ] should also be taken into account for the dust particles . in these tests , the more restrictive conditions will occur at the shock front . to set up the shock tube test , an ensemble of particles with @xmath205 , @xmath206 , @xmath207 , @xmath208 , are evenly distributed in a one - dimensional bounded domain @xmath168 . to model the density jump , a different number of particles is used at each side of the discontinuity . in particular , since it is a one - dimensional case @xmath209 where @xmath144 is the total number of particles . particle masses are calculated as in the previous sections . [ fig : test3 ] , presents the results for two weakly dragged cases in a constant drag regime . the left panel represents a case with @xmath210 , while the right panel represents non linear drag regime ( equations [ kst ] and [ cd3 ] ) with @xmath211 . despite not having an analytical solution for the transient phase , in both figures , the long - term analytical solution of the problem has been added ( dotted line ) as a guideline . in both cases , the obtained solution compares very favourably with the results previously obtained by @xcite through the use of explicit / implicit methods . in the left panel of fig . [ fig : test3b ] , a strongly dragged case with @xmath212 is presented for two different resolutions in the non - linear regime . in this case , the analytical solution is known ( solid line ) , and as can be seen , it is well matched by the numerical results if enough resolution is used . in the right panel of fig .  [ fig : test3b ] the same case is presented for a @xmath157 case with 569 gas and 50 dust particles . as can be seen , despite the reduced dust resolution , the correct result is obtained and there is no evidence of overdissipation .    . such a calculation would be impossible with an explicit time step.,width=302 ]    it is also very interesting to compare the results obtained using the semi - implicit method , with those obtained with an explicit integration scheme , for a very high drag regime with @xmath213 . as can be seen in fig . [ fig : test3c ] , whereas excess dissipation in the explicit calculation gives an incorrect solution , our semi - implicit method avoids the problem . as previously mentioned , the source of overdissipation in the semi - implicit method comes from the incapacity of the algorithm to estimate the local barycentric velocity , due to the lack of resolution . however , and in contrast with an ordinary explicit method , if the semi - implicit method is used , the error in the barycentric velocity estimation is only committed once per gas integration time - step . on the contrary , if an explicit method is used , due to the courant condition of the drag interaction , the error in the estimation of the drag acceleration is committed a lot more times per gas integration time - step , leading to a very poor result . this effect will also occur if a conventional implicit integration scheme is used ( e.g. * ? ? ? * ) . finally , in the uppermost panels of fig . [ fig : test4a ] the result of a sedov blast test with @xmath184 , @xmath179 , and @xmath214 ( @xmath215 ) is presented . in this case , the gas integration time - step is set by the courant time condition at the shock front . in the test , a total of @xmath216 particles are evenly distributed in a three - dimensional cartesian grid with @xmath217 . particle masses are calculated as in the previous sections . the dust grid is displaced with respect to the gas one by half of the gas particle separation in each direction . to model the explosion a total thermal energy of @xmath218 code units is distributed over the particles inside a certain radius ( @xmath219 ) . as a comparison , the sedov blast test performed with an ordinary explicit integration scheme , is also presented in the lower panels of fig . [ fig : test4a ] . as can be seen , there are no significant differences . .computational time increase factors as a function of the drag strength in the sedov test . the computational time of each simulation with the explicit method is divided by the computational time of the semi - implicit method . in the semi - implicit method , since the time - step of the simulation is exclusively determined by the gas courant time condition , and the dust - to - gas ratio is small , no noticeable extra computational effort is needed if the drag strength is increased . [ cols="^,^",options="header " , ]     again , in this case , no evidence of the resolution limitation has been found , due to the lower density ratio between the gas and dust components . since the typical gas - to - dust ratios in the interstellar medium are very similar to the ones used in the sedov test , we expect the method to be useful in realistic astrophysical simulations . additionally , we have used this test to compare the computational time of the method , with that of a traditional explicit integration . in table 1 , a comparison of the computational time for several different drag strengths , is presented for both cases . in each case , the computational time spent by each simulation is divided by the computational time of the semi - implicit method in the @xmath215 case . as can be seen , as the drag strength is increased , the explicit integrator computational time is increased , by several orders of magnitude , with respect to the computational time spent by semi - implicit method for @xmath215 . on the contrary , the computational time of the semi - implicit method remains stable , since the integration time - step is exclusively determined by the gas courant condition , and is independent of the drag strength . as can be seen in fig . [ fig : test4b ] , arbitrarily large values for the drag coefficient can be used . this value would be completely prohibitive in any two - fluid explicit integration method . the final test performed was the mid - plane settling of dust particles in a one - dimensional vertical section of an isothermal disk with @xmath169 and @xmath170 . to set up the test , 100 gas and 100 dust particles , with @xmath184 and @xmath220 , are evenly distributed over a one - dimensional domain ( @xmath221 ) . particle masses are assigned following the same procedure as in the previous sections . an external acceleration @xmath222 is used to simulate the vertical component of the gravitational field from the star at the centre of the disk , where @xmath223 is the angular frequency . ( see appendices a and b for a detailed explanation about how to implement external forces in the integration scheme ) . no boundaries have been used , and since one does not expect shocks to be important , the use of artificial viscosity is avoided . the evolution equations of the system are given by @xmath224 @xmath225 in order for the system to relax , the gas particles are evolved under gravitational and pressure forces , until the hydrostatic equilibrium condition is attained . whenever hydrostatic equilibrium is reached , equations [ eu_disk1 ] and [ eu_disk2 ] can be solved @xmath226 @xmath227 giving the gas hydrostatic density profile : @xmath228 which is valid as long as @xmath229 . in fig . [ fig : test5a ] , the initial gas density profile of the isothermal disk is presented . as can be seen , the gas perfectly reproduces a gaussian density profile with @xmath230 , @xmath231 and @xmath232 . after gas relaxation , drag forces are switched on , and evolution is started again . if the drag coefficient @xmath233 is high enough , dust particles reach a limiting velocity , given by the solution of equations [ eq1 ] and [ eq2 ] . @xmath234    in fig . [ fig : test5b ] , the dust component velocity as a function of @xmath235 is presented for two cases ( @xmath236 and @xmath237 ) . as can be seen , the correct limiting velocity of the dust component is reached in both cases . because @xmath238 , the momentum transferred between the dust and gas phases is rather small , and the gas component remains very close to the hydrostatic equilibrium . as can be seen in the @xmath237 case ( right plot ) , dust particles almost instantaneously reach its limiting velocity . on the contrary , if @xmath236 ( left plot ) , particles need more time to reach the limiting velocity and the transitory state can be seen for @xmath239 . in order to check whether the algorithm is capable of correctly reproducing such a transitory regime , the velocity as a function of @xmath235 for a single sph dust particle can be compared with the numerical solution of equations [ eu_disk1 ] and [ eu_disk2 ] . in fig . [ fig : test5c ] the evolution of a single sph dust particle is plotted for three different @xmath240 values . circles represent the velocity of the particle , for different time steps , as it falls down towards the disk mid - plane . dashed lines represent the numerical solution of equations [ eu_disk1 ] and [ eu_disk2 ] for each case , while solid lines represent the limiting velocity for each case as given by equation [ v_lim ] . if @xmath241 , the particle does not have time to reach the limiting velocity , and simply suffers velocity damping while it oscillates around the disk mid - plane . as can be seen , a perfect agreement is achieved with the theoretical behaviour . if @xmath236 , the dust particle reaches the limiting velocity at @xmath242 , in perfect agreement with the numerical solution of equations [ eu_disk1 ] and [ eu_disk2 ] , and explaining the global velocity profile of the dust component ( as seen in fig . [ fig : test5b ] ) . for @xmath237 , although the theoretical solution is approximately obtained , some oscillations of the particle velocity can be observed . such oscillations occur due to the low number of gas particles present in the outermost parts of the disk . if a higher resolution simulation is performed ( 1000 gas particles ) , the oscillations disappear , and the velocity of the dust particle closely matches the analytical solution . .,width=302 ]          + a new method has been proposed to avoid explicit integration of the time evolution equations of small dust grains in the two fluid sph approach . through the use of semi - analytic solutions for the decay of the gas and dust relative velocity , the present method has been able to reproduce all the features of the previous two fluid sph approach of @xcite , with the advantage of a considerable gain in computational time in strong drag regimes . due to its strictly dissipative nature , the velocity changes induced by the drag force can be estimated without the need for explicit acceleration recalculations or iterative procedures , even when the stopping time becomes much shorter than the gas evolutionary time - scale . the method is numerically stable , and always provides convergence towards the analytical solutions as the resolution is increased . the method has also been capable of reproducing the correct behaviour of the drag force for all regimes . in the weak drag regime , the method is theoretically equivalent to a standard explicit integration , both in accuracy and computational efficiency , as long as strong gradients are not present in the immediate neighbourhood of dust particles . in the high drag regime , the method is capable of reproducing all the expected features of dust / gas mixtures . the results obtained in the test cases are completely analogous to those found by @xcite through the use of standard explicit and implicit methods .    in agreement with previous studies @xcite , a resolution limit has been found for the method in the dustywave experiment . for high drag regimes with dust - to - gas ratios of order unity , the resolution should exceed @xmath195 in order to avoid overdissipation . however , in the shock tube experiment , our method avoids the effects of overdissipation , which until now has been considered to be one of the main limitations of the two - fluid sph approach . furthermore , it has also been demonstrated that for low drag regimes , and even for high drag regimes with low dust - to - gas ratios , the number of dust particles present in the simulation becomes irrelevant , and the accuracy of the solution is only dependent on having sufficient gas resolution . since in the vast majority of astrophysical applications the dust - to - gas ratio is expected to be rather low , only a good gas resolution will be necessary to avoid overdissipation . however , special attention must be payed to this limitation , since it will be very difficult to completely avoid overdissipation in complex global simulations , especially if one expects abrupt changes in the dust - to - gas ratios . we thank the anonymous referee , whose thoughtful and thorough report not only resulted in substantial improvements to the paper , but as a by - product also meant that we improved the numerical algorithm itself . we also thank joe monaghan and daniel price for very useful discussions . figures 2 to 9 have been created using splash @xcite , a sph visualization tool publicly available at http://users.monash.edu.au/@xmath194dprice/splash . the calculations for this paper were performed on the dirac complexity machine , jointly funded by stfc and the large facilities capital fund of bis , and the university of exeter supercomputer , a dirac facility jointly funded by stfc , the large facilities capital fund of bis and the university of exeter . this work was also supported by the stfc consolidated grant st / j001627/1 . 23    armitage p.  j.,2010 , cambridge university press    ayliffe b.  a. , laibe g. , price d.  j. , bate m.  r , 2011 , mnras , 423 , 1450    bate m. , 1995 , phd thesis , univ . cambridge    fehlberg e. , low - order classical runge - kutta formulas with step size control and their application to some heat transfer problems , nasa technical report 315    fulk d.a . , quinn d.w . , 1996 , j. comput . phys . , 17 , 19    laibe g. , price d.  j. , 2011 , mnras , 418 , 1491    laibe g. , price d.  j. , 2012a , mnras , 420 , 2345    laibe g. , price d.  j. , 2012b , mnras , 420 , 2365    laibe g. , price d.  j. , 2014 , mnras , 440 , 2136    laibe g. , price d.  j. , 2014 , mnras , 440 , 2147    marble , f. , 1970 , annual review of fluid mechanics , 2 , 397    monaghan j. j. , 1992 , annual review of astronomy and astrophysics , 30 , 543    monaghan j. j. , 1997 , journal of computational physics , 136 , 298    monaghan j. j. , 2002 , mnras , 335 , 843    monaghan j. j. , 1997 , journal of computational physics , 138 , 801    monaghan j. j. , kocharyan a. , 1995 , computer physics communications , 87 , 225    price d.  j. , 2007 , publications of the astronomical society of australia , 24 , 159    price d.  j. , 2008 , journal of computational physics , 227 , 10040    price d.j . , monaghan j. j. , 2004 , mnras , 348 , 139    randles p.w . , and libersky l.d . , 1996 , smoothed particle hydrodynamics some recent improvements and applications , computer methods in applied mechanics and engineering , 138 , 375    saffman p.g.,on the stability of laminar flow of a dusty gas , j. fluid mechanics , 13 , 120    serna a. , alimi j .- m . , chieze j .- , 1995 , astrophysical journal , 461 , 884    sedov l. , 1959 , similarity and dimensional methods in mechanics , new york : academic    springel v. , hernquist l. , 2002 , mnras , 333 , 649    weidenschilling , s.j . , 1977 , mnras , 180 , 57    whipple , f.l . , 1972 , from plasma to planet , ed a. elvius , wiley , london , p.211    wtezstein m. , nelson a.  f. , naab t. , burkert a. , 2009 , apjs , 184 , 298 in order to capture both the gas and dust evolution described by equations [ eu1 ] , [ eu2 ] , [ eu3 ] , [ cont1 ] , and [ cont2 ] , the dust - gas drag equations must be coupled with an explicit hydrodynamical integrator . combination with an integrator incorporating the gas pressure gradients is needed . one chosen integrator is a second order runge - kutta - fehlberg scheme and the combined scheme can be summarized as follows @xmath243 for the first half of the time - step , and @xmath244 for the full time - step . we have introduced the dust and gas components external accelerations @xmath245 and @xmath246 , in order to account for forces like gravity , physical viscosity or radiation pressure . this method relies on the possibility of considering pressure and drag forces as separable interactions . as the performed tests have shown , it seems to be a good assumption . + another particularly useful property of the present method is its capacity to predict the correct modified sound speed of the dust / gas mixture , as a function of the dust / gas ratio . by substituting the pre - dragged quantities @xmath247 and @xmath248 and the expression for the @xmath101 parameter into the @xmath249 and @xmath250 equations , one can convert the two - step method into an equivalent one - step method given by the set of equations for the first half of the time - step @xmath251 where @xmath252 , @xmath253 and @xmath254 for the full time - step , where we have defined @xmath255 as we can see , in this set of equations , dust can be no longer considered pressureless . it suffers an acceleration due to pressure gradient , and possesses an effective density @xmath256 . this result can be understood if one realizes that a purely dissipative force does not always lead to a velocity decrease . because drag is a purely dissipative force , it will always lead to a decrease in the relative velocity between dust and gas components . but sometimes , the only way to decrease such a relative velocity is to accelerate the dust component . also , the effective densities @xmath256 and @xmath257 can be understood as the effective inertial response of the dust and gas components to the effective pressure terms . the weaker the drag force is , the higher the pressure gradient must be to accelerate the dust component . if @xmath258 , @xmath259 , @xmath260 , and the equations for the change in velocity of the dust and gas components become @xmath261 for the first half time - step and @xmath262 for the full time - step . the effective dust density term @xmath256 has become infinitely big , so the dust does not respond at all to the pressure gradient terms . that is , gas and dust decouple , and gas evolves as a single component fluid with sound speed @xmath263 . if , on the contrary , @xmath264 , @xmath265 , @xmath266 , and the equations for the change in velocity of the dust and gas components become this time @xmath267 for the first half time - step and @xmath268 for the full time - step . both effective density terms @xmath256 and @xmath257 become equal , so both dust and gas components evolve as a single component fluid , with the total mass of the mixture being advected however , and since only the gas component can produce real pressure , they travel with a modified sound speed @xmath269 , exactly as predicted by theory ( see for example @xcite ) . as can be seen in equations a9 , both phases adopt in this regime the barycentric velocity in just one time - step , as it corresponds to a case where @xmath270 . despite being particularly useful to visualize the behaviour of dust and gas mixtures , and to show that the effective sound speed of the mixture is the expected one in the strong drag regime , we still recommend using the two - step method given by equations a1 and a2 . it is clearly technically easier to implement into a pre - existing sph code . the second chosen integrator was a modification of the predictor - corrector scheme of @xcite and it can be summarized as follows @xmath271 for the predictor phase and @xmath272\delta t/2,\\ \tilde{\textbf{v}}^{t+1}_{\rm g } & = \textbf{v}^{t+1/2}_{\rm g } - \left[\left(\frac{\mathbf{\nabla } p}{\rho_{\rm g}}\right)_{t+1/2 } - \left(\frac{\mathbf{\nabla } p}{\rho_{\rm g}}\right)_{t}\right]\delta t/2 + \left[\textbf{a}_{\rm ext , g}^{t+1/2 } -   \textbf{a}_{\rm ext , g}^{t}\right]\delta t/2 , \\ \textbf{v}^{t+1}_{\rm d } & = \tilde{\textbf{v}}^{t+1}_{\rm d } -\xi\left[(\tilde{\textbf{v}}^{t+1}_{\rm d}-\tilde{\textbf{v}}^{t+1}_{\rm g})-(\textbf{v}^{t+1/2}_{\rm d}-\textbf{v}^{t+1/2}_{\rm g})\right],\\ \textbf{v}^{t+1}_{\rm g } & = \tilde{\textbf{v}}^{t+1}_{\rm g } + \frac{\rho_{\rm d}}{\rho_{\rm g}}\xi\left[(\tilde{\textbf{v}}^{t+1}_{\rm d}-\tilde{\textbf{v}}^{t+1}_{\rm g})-(\textbf{v}^{t+1/2}_{\rm d}-\textbf{v}^{t+1/2}_{\rm g})\right],\\ \textbf{r}^{t+1}_{\rm d } & = \textbf{r}^{t+1/2}_{\rm d } + ( \textbf{v}^{t+1}_{\rm d}-\textbf{v}^{t+1/2}_{\rm d})\delta t/3,\\ \textbf{r}^{t+1}_{\rm g } & = \textbf{r}^{t+1/2}_{\rm g } + ( \textbf{v}^{t+1}_{\rm g}-\textbf{v}^{t+1/2}_{\rm g})\delta t/3 , \end{aligned}\ ] ] for the corrector phase .
a method to avoid the explicit time integration of small dust grains in the two fluid gas /dust smoothed particle hydrodynamics ( sph ) approach is proposed . by assuming a very simple exponential decay model for the relative velocity between the gas and dust components , all the effective characteristics of the drag force can be reproduced . 
 a series of tests has been performed to compare the accuracy of the method with analytical and explicit integration results . 
 we find that the method performs well on a wide range of tests , and can provide large speed ups over explicit integration when the dust stopping time is small . 
 we have also found that the method is much less dissipative than conventional explicit or implicit two - fluid sph approaches when modelling dusty shocks . 
 [ firstpage ]    hydrodynamics - methods : numerical - planets and satellites : formation - protoplanetary discs - dust , extinction .