# An implicit algorithm for simulating the dynamics of small dust grains with
smoothed particle hydrodynamics

Daniel Elsender and Matthew R. Bate
Department of Physics and Astronomy, University of Exeter, Stocker Road,
Exeter EX4 4QL, UK E-mail<EMAIL_ADDRESS>(DE)

(Accepted XXX. Received YYY; in original form ZZZ)

###### Abstract

We present an implicit method for solving the diffusion equation for the
evolution of the dust fraction in the terminal velocity approximation using
dust-as-mixture smoothed particle hydrodynamics (SPH). The numerical scheme
involves casting the dust diffusion equation into implicit form, rearranging
into its resolvent cubic equation and solving analytically. This method is
relevant for small grains that are tightly coupled to the gas, such as sub-
micron dust grains in the interstellar medium or millimetre-sized dust grains
in protoplanetary discs. The method avoids problems with the variable used to
evolve the dust fraction becoming negative when evolved explicitly and is fast
and accurate, avoiding the need for dust stopping time limiters and
significantly reducing computational expense. Whilst this method is an
improvement over using the explicit terminal velocity approximation method, as
with any dust-as-mixture method it still fails to give accurate solutions in
the limit of large (weakly coupled) grains.

###### keywords:

hydrodynamics – methods: numerical – protoplanetary discs – dust, extinction –
ISM: kinematics and dynamics

††pubyear: 2023††pagerange: An implicit algorithm for simulating the dynamics
of small dust grains with smoothed particle hydrodynamics–A

##  1 Introduction

Dust emission is the dominant signature in cold astrophysical observations.
Small grains absorb ultraviolet/optical light emitted from stars and re-emit
this energy in the infrared. Instruments such as the mid-infrared instrument
(MIRI) onboard the James Webb Space Telescope (JWST) is primed to detect these
emissions from dust. Infrared observations can be used in conjunction with
radio observations, such as those using the Atacama Large (sub)Millimeter
Array (ALMA), to explore the structures of protoplanetary discs. Within
protoplanetary discs it is dust that provides, at least in part, the material
from which planets form. In addition to this, dust grains collide with gas
particles, altering the dynamical evolution of both phases. Understanding the
dynamics of dust in these environments is crucial for the interpretation of
observations that rely upon assumptions about the dust-to-gas ratio, as dust
structures can be vastly different to structures in the gas.

Smoothed particle hydrodynamics (SPH) is a popular numerical method for
modelling astrophysical fluid dynamics. Previously dust and gas mixtures have
been modelled in SPH using separate particles for the dust and gas, originally
by Monaghan & Kocharyan (1995), where the dust is coupled to the gas via drag.
Such methods are sometimes referred to as two-fluid methods, but are less
ambiguously referred to as dust-as-particles. This dust-as-particles method
was applied in the context of protoplanetary discs by Maddison et al. (2003)
and Barrière-Fouchet et al. (2005). Subsequently, Laibe & Price (2012a, b)
presented a new method for simulating dust-as-particles with SPH. Grid codes
have also been adapted to model dust-as-particles. Bai & Stone (2010)
introduced a hybrid particle-gas scheme to simulate dust-as-particles in dusty
gas mixtures with the grid code athena; they also present a fully implicit
solver. The grid-code fargo3d uses an implicit numerical scheme to model
multiple species of dust (Benítez-Llambay et al., 2019). Implicit timestepping
methods have also been employed in SPH dust-as-particle methods. Monaghan
(1997) modified this method to treat the drag terms implicitly via pairwise
interactions. More recently, Laibe & Price (2012b) introduced an implicit
timestepping method for their dust-as-particle method based upon the pairwise
interaction method of Monaghan (1997). Lorén-Aguilar & Bate (2014, 2015)
developed a semi-implicit approach to dust-as-particle methods in SPH,
although this scheme is limited to the use of global timesteps. Stoyanovskaya
et al. (2018); Stoyanovskaya et al. (2020) presented an implicit dust-as-
particles SPH method for both linear and non-linear drag regimes.

Modelling the dust as a set of particles in essence doubles the memory
required to perform the calculation and usually significantly increases the
time to completion compared to an equivalent calculation with only gas. In the
past decade much effort has gone into understanding and modelling the dynamics
of small dust grains, that are highly coupled with the gas phase. In a series
of papers Laibe & Price (2014a, b, c) discuss the limitations of the dust-as-
particles method in the high drag (small grain) regime. To circumvent these
issues they developed a dust-as-mixture method (sometimes called one-fluid)
for modelling highly coupled dusty gas mixture by framing the equations in the
barycentric reference frame of the fluid. Price & Laibe (2015) presented a
simplified version by adopting the terminal velocity approximation (Youdin &
Goodman, 2005) and modelling the dust fraction evolution as a diffusion-like
equation. The dust diffusion method is far less computationally expensive than
an explicit dust-as-particles method in the limit of small grains. In addition
to being deployed in SPH, the dust-as-mixture method for modelling small grain
dynamics has been implemented in the adaptive mesh refinement code RAMSES
(Lebreuilly et al., 2019).

The dust-as-mixture method of Price & Laibe (2015) has been extended to model
the dynamics of multiple small dust grain populations by Hutchison, Price &
Laibe (2018). This dust-as-mixture method is an effective way to simulate the
dynamics of small dust grains in the limit of the terminal velocity
approximation. However, although the method enforces positivity on the dust
fraction, $\epsilon$, (by evolving a variable $s=\sqrt{\epsilon\rho}$), where
$\rho$ is the total mass density) it does not limit the dust fraction to less
than unity. This issue was solved by Ballabio et al. (2018) who use a
different parameterisation of dust fraction in their evolution equations.
Namely the dust fraction is defined as

$\epsilon=\frac{s^{2}}{1+s^{2}},$ (1)

and it is the dust variable $s$ that is evolved. An alternative formulation
using a dust variable defined by $\epsilon=\sin^{2}(s)$ was used by Hutchison
et al. (2022). This formulation also enforces $\epsilon\in[0,1]$ and was found
to produce more accurate solutions than equation 1, but it is more
computationally expensive due to the need to evaluate trigonometric functions.

For all of the above formulations, while the positivity of the dust fraction
is guaranteed, there is nothing stopping the dust variable, $s$, itself from
becoming negative. In the regions of simulations where there are steep dust
fraction gradients, while enforcing exact dust mass conversation, the dust
variable can undershoot zero and become negative. If it does so this can cause
the dust fraction to become positive in these regions where, otherwise, there
should be no dust. This consequence of the mathematics of the model can be
thought of as dust leaking unphysically into these dust sparse regions. This
can become a significant issue when modelling dust grains that are moderately
coupled to the gas, for example, in protoplanetary discs where larger dust
grains would be expected to settle to the midplane of the disc they can ‘leak’
into or be retained in the upper, and lower, regions of the disc relative to
the midplane. This has been avoided in the past by sacrificing exact dust mass
conservation and simply neglecting any change in the dust fraction or dust
variable that would make it negative (e.g., see section 4.4.2 of Price &
Laibe, 2015). The resulting dust mass conservation is better when using
equation 1 or the $\epsilon=\sin^{2}(s)$ formulations than when using
$s=\sqrt{\epsilon\rho}$ or evolving the dust fraction directly, but it is
still not ideal.

An additional limitation of the dust-as-mixture method is that it requires a
timestep criterion whereby timesteps become small in regions of low dust
fraction, or where dust grains become less well coupled (i.e., large). To
prevent calculations having prohibitively short timesteps, Ballabio et al.
(2018) introduced a mass flux limiter to stop rapid dust diffusion for
particles in regions of high dust fraction gradients by artificially limiting
the stopping time. The problem with this is that it unphysically modifies the
drag acting on the dust during a simulation and, although it provides a
numerical solution, the degree to which the solution resembles reality is
unclear.

In this paper, we present an implicit method for solving the dust diffusion
equation of Price & Laibe (2015), using the dust variable defined by equation
1 (Ballabio et al., 2018). This avoids the need for a dust timestep criterion
and negates the need for a stopping time limiter. In addition, it better
handles the evolution of the dust variable in regions of large dust gradients
due to the iterative nature of the algorithm. Depending on the dust grain
sizes, the method also significantly speeds up the time to completion of
three-dimensional simulations of protostellar collapse and protoplanetary disc
evolution simulations, and yet yields equivalent results to those obtained
using the fully explicit method.

Our aim with this work is to improve calculations in the limit of small,
coupled grains. Although we have applied this to the particle-as-mixture
method that uses the terminal velocity approximation, a similar method could
be derived to evolve the dust-as-mixture equations without the additional
terminal velocity approximation. Applying such an implicit method to the full
dust-as-mixture equations of Laibe & Price (2014a, b) would in principle allow
grains with long stopping times to be more accurately modelled. This may be
worth investigating in the future, but it is beyond the scope of this paper.
Furthermore, while the full set of dust-as-mixture equations can provide a
better model for intermediate-sized grains than applying the terminal velocity
approximation, the fundamental assumption of a particle-as-mixture method that
the velocity field is single valued does eventually break down, as has been
shown by Laibe & Price (2014b) and Bate & Loren-Aguilar (2017). Finally, we
note that Hutchison et al. (2016) used an implicit method to evolve the the
differential velocity between gas and dust in full dust-as-mixture SPH method,
but they still solved the dust diffusion equation explicitly.

##  2 Method

###  2.1 The dust-as-mixture dust method

In the derivation of the dust-as-mixture formulation (Laibe & Price, 2014a),
the continuum equations are rewritten in the barycentric reference frame of
the fluid. This involves replacing the velocities of the gas and dust phases
with a single barycentric velocity of the mixture,

$\boldsymbol{v}=\frac{\rho_{\mathrm{g}}\boldsymbol{v}_{\mathrm{g}}+\rho_{\mathrm{d}}\boldsymbol{v}_{\mathrm{d}}}{\rho_{\mathrm{g}}+\rho_{\mathrm{d}}},$
(2)

where $\boldsymbol{v}$ is the velocity, $\rho$ is the density and the
subscripts d and g denote the dust and gas quantities, respectively. The
differential velocity between the two phases is defined as,

$\Delta\boldsymbol{v}=\boldsymbol{v}_{\mathrm{d}}-\boldsymbol{v}_{\mathrm{g}}.$
(3)

The total density of the mixture is,

$\rho=\rho_{\mathrm{g}}+\rho_{\mathrm{d}},$ (4)

and the dust fraction is defined as,

$\epsilon\equiv\rho_{\mathrm{d}}/\rho.$ (5)

The dust-as-mixture equations can then be written in the form

$\displaystyle\derivative{\rho}{t}=-\rho\left(\divergence{\boldsymbol{v}}\right),$
(6)
$\displaystyle\derivative{\epsilon}{t}=-\frac{1}{\rho}\divergence{\left[\epsilon\left(1-\epsilon\right)\rho\Delta\boldsymbol{v}\right]},$
(7)
$\displaystyle\derivative{\boldsymbol{v}}{t}=-\frac{\gradient{P}}{\rho}-\frac{1}{\rho}\divergence{\left[\epsilon\left(1-\epsilon\right)\rho\Delta\boldsymbol{v}^{2}\right]}+\boldsymbol{f},$
(8)
$\displaystyle\derivative{\Delta\boldsymbol{v}}{t}=-\frac{\Delta\boldsymbol{v}}{t_{\mathrm{s}}}+\frac{\gradient{P}}{\left(1-\epsilon\right)\rho}-\left(\Delta\boldsymbol{v}\cdot\gradient{}\right)\boldsymbol{v}+\frac{1}{2}\divergence{\left[\left(2\epsilon-1\right)\Delta\boldsymbol{v}\Delta\boldsymbol{v}\right]},$
(9)
$\displaystyle\derivative{u}{t}=-\frac{P}{\rho_{\mathrm{g}}}\divergence{\left(\boldsymbol{v}-\epsilon\Delta\boldsymbol{v}\right)}+\epsilon\Delta\boldsymbol{v}\cdot\nabla
u+\epsilon\frac{\Delta\boldsymbol{v}^{2}}{t_{\mathrm{s}}},$ (10)

where $P$ is the gas pressure, $\boldsymbol{f}$ represents non-fluid forces
acting on the mixture (e.g., gravity), $u$ is the internal energy of the gas,
and $t_{\mathrm{s}}$ is the stopping time of the dust.

The stopping time when the dust is in the Epstein drag regime can be expressed
as

$t_{\mathrm{s}}=\frac{\hat{\rho}_{\mathrm{s}}r_{\mathrm{s}}}{\rho_{\mathrm{g}}v_{\mathrm{th}}},$
(11)

where $\hat{\rho}_{\mathrm{s}}$ is the intrinsic grain density,
$r_{\mathrm{s}}$ is the grain radius, and the velocity of the gas molecules
due to thermal motion is

$v_{\mathrm{th}}=\sqrt{\frac{8k_{\mathrm{B}}T}{\pi\mu m_{\mathrm{H}}}},$ (12)

where $k_{\mathrm{B}}$ is Boltzmann’s constant, $T$ is the gas temperature,
$\mu$ is the mean molecular weight of the gas, and $m_{\mathrm{H}}$ is the
atomic mass of hydrogen.

This formalism can be simplified when the stopping time is small compared to
the hydrodynamic time scale (Price & Laibe, 2015). In SPH terms this means
$t_{\mathrm{s}}<h/c_{\mathrm{s}}$, where $h$ is the SPH particle smoothing
(resolution) length and $c_{\mathrm{s}}$ is the sound speed. When this is
satisfied, we are in the terminal velocity regime (Youdin & Goodman, 2005),
i.e. the relative velocity of the dust and gas reaches terminal velocity
quickly due to the drag and pressure forces balancing. In this regime the
dust-as-mixture equations can be simplified by ignoring the time dependence of
the differential velocity, and any terms that are second order in
$t_{\mathrm{s}}$. Then

$\Delta v=t_{\mathrm{s}}\frac{\gradient{P}}{\rho_{\mathrm{g}}}.$ (13)

The continuum equations then become

$\displaystyle\derivative{\rho}{t}=-\rho\left(\divergence{\boldsymbol{v}}\right),$
(14)
$\displaystyle\derivative{\epsilon}{t}=-\frac{1}{\rho}\divergence{\left(\epsilon
t_{\mathrm{s}}\gradient{P}\right)},$ (15)
$\displaystyle\derivative{v}{t}=-\frac{\gradient{P}}{\rho}+\boldsymbol{f},$
(16)
$\displaystyle\derivative{u}{t}=-\frac{P}{\rho_{\mathrm{g}}}\left(\divergence{\boldsymbol{v}}\right)+\frac{\epsilon
t_{\mathrm{s}}}{\rho_{\mathrm{g}}}\left(\gradient{P}\cdot\gradient{u}\right).$
(17)

The inclusion of equation 15 means an additional time step criterion is
required. Price & Laibe (2015) provide such a criterion in the form

$\Delta t<C_{\mathrm{0}}\frac{h^{2}}{\epsilon
c_{\mathrm{s}}^{2}t_{\mathrm{s}}},$ (18)

where $C_{\mathrm{0}}$ is the Courant number. This implies that the time steps
are constrained when the stopping times are long, i.e. when drag is low / dust
grains are not well coupled to the gas. These equations can be discretised
into the SPH form. Notably, the continuity and momentum equations have the
same discretisation as the gas-only SPH equations. We must have new
discretisations for the dust fraction, and the thermal energy evolution
equations.

To discretise the dust fraction evolution equation whilst maintaining
positivity of the dust fraction, and limiting the fraction to below unity,
thus prevents the dust fraction from becoming unphysical, we use the
parameterisation of equation 1 introduced by Ballabio et al. (2018) resulting
in the dust evolution variable

$s=\sqrt{\frac{\epsilon}{1-\epsilon}}.$ (19)

We compute the time evolution of this variable

$\derivative{s}{t}=\frac{1}{2s\left(1-\epsilon\right)^{2}}\derivative{\epsilon}{t},$
(20)

using a direct second derivative (Price & Laibe, 2015). Substituting Eq. 15
into Eq. 20 yields

$\derivative{s}{t}=-\frac{1}{2\rho(1-\epsilon)^{2}}\\{\divergence{[s(1-\epsilon)t_{\mathrm{s}}\gradient{P}]}+(1-\epsilon)t_{\mathrm{s}}\gradient{P}\cdot\gradient{s}\\}.$
(21)

Using the implementation of Ballabio et al. (2018), the SPH discretisation of
this equation is

$\displaystyle\derivative{s_{i}}{t}=-\frac{1}{2\rho_{i}(1-\epsilon_{i})^{2}}\sum_{j}\left[\frac{m_{j}s_{j}}{\rho_{j}}(D_{i}+D_{j})(P_{i}-P_{j})\frac{\bar{F}_{ij}}{r_{ij}}\right],$
(22)

where $D_{i}=t_{s,i}(1-\epsilon_{i})$,
$\bar{F}_{ij}=\left[F_{ij}(h_{i})+F_{ij}(h_{j})\right]/2$,
$r_{ij}=|\boldsymbol{r}_{ij}|$, and
$\mathbf{\hat{r}}_{ij}F_{ij}(h_{i})=\gradient_{i}W_{ij}(h_{i})$, where
$W_{ij}(h_{i})$ is the SPH kernel. This work utilises the M6 quintic spline
kernel.

The energy of the system can be written as

$E=\sum_{i}m_{i}\left[\frac{1}{2}v^{2}_{i}+(1-\epsilon_{i})u_{i}\right],$ (23)

which is conserved if

$\derivative{E}{t}=\sum_{i}m_{i}\left[v_{i}\derivative{v_{i}}{t}+(1-\epsilon_{i})\derivative{u_{i}}{t}-u_{i}\derivative{\epsilon_{i}}{t}\right]=0.$
(24)

Removing the non-dust terms, assuming they conserve energy, and substituting
in our parameterisation for dust fraction gives

$\displaystyle\evaluated{\derivative{E}{t}}_{\text{missing}}{dust}$
$\displaystyle=\sum_{i}m_{i}(1-\epsilon_{i})\evaluated{\derivative{u_{i}}{t}}_{\text{missing}}{dust}$
$\displaystyle\quad\quad-\sum_{i}m_{i}u_{i}\frac{2s_{i}}{\left(1+s_{i}^{2}\right)^{2}}\derivative{s_{i}}{t}=0,$
(25)

using Eq. 22 and rearranging yields

$\begin{split}\evaluated{\derivative{u_{i}}{t}}_{\text{missing}}{dust}=-\frac{1}{2(1-\epsilon_{i})}\sum_{j}&m_{j}\frac{s_{i}s_{j}}{\rho_{j}}\left(D_{i}+D_{j}\right)\\\
&\left(P_{i}-P_{j}\right)\left(u_{i}-u_{j}\right)\frac{\bar{F}_{ij}}{r_{ij}}.\end{split}$
(26)

When incorporated into the usual SPH equations for a gas, equation 14 can be
expressed as

$\derivative{\rho_{i}}{t}=\sum_{j}m_{j}(\mathbf{v}_{i}-\mathbf{v}_{j})\cdot\gradient_{i}W_{ij}(h_{i}),$
(27)

although this equation is not explicitly solved as the density is computed by
the usual weighted sum over neighbours, and equations 16, and 17 become

$\displaystyle\derivative{\mathbf{v}_{i}}{t}$
$\displaystyle=-\sum_{j}m_{j}\left(\frac{P_{i}}{\rho_{i}^{2}\Omega_{i}}\gradient_{i}{W_{ij}(h_{i})}+\frac{P_{j}}{\rho_{j}^{2}\Omega_{j}}\gradient_{j}{W_{ij}}(h_{j})\right)$
$\displaystyle+\boldsymbol{f}+\Pi_{\text{AV}},$ (28)
$\displaystyle\derivative{u_{i}}{t}$
$\displaystyle=\frac{P_{i}}{\Omega_{i}\rho_{i}^{2}}\sum_{j}m_{j}(\mathbf{v}_{i}-\mathbf{v}_{j})\cdot\gradient_{i}{W_{ij}(h_{i})}+\left(\derivative{u_{i}}{t}\right)^{\text{diss}}$
$\displaystyle-\frac{1}{2(1-\epsilon_{i})}\sum_{j}m_{j}\frac{s_{i}s_{j}}{\rho_{i}\rho_{j}}\left(D_{i}+D_{j}\right)\left(P_{i}-P_{j}\right)\left(u_{i}-u_{j}\right)\frac{\bar{F}_{ij}}{r_{ij}},$
(29)

where $\Omega_{i}$ is the term related to the smoothing length gradients
defined as

$\Omega_{i}\equiv
1-\partialderivative{h_{i}}{\rho_{i}}\sum_{j}m_{j}\partialderivative{W_{ij}(h_{i})}{h_{i}}.$
(30)

The $\Pi_{AV}$ term provides shock capturing via artificial viscosity, which
in the SPH code we used is formulated as

$\Pi_{ij}=-\sum_{j}\frac{m_{j}}{\rho_{ij}}\left[\frac{q_{ij}^{i}}{\rho_{i}\Omega_{i}}\gradient_{i}{W_{ij}(h_{i})}+\frac{q_{ij}^{j}}{\rho_{j}\Omega_{j}}\gradient_{j}{W_{ij}(h_{j})}\right],$
(31)

where $q_{ij}$ is defined as

$q^{i}_{ij}=\begin{cases}-\frac{1}{2}\rho_{i}v_{\text{sig},i}\mathbf{v}_{ij}\cdot\mathbf{\hat{r}}_{ij},&\mathbf{v}_{ij}\cdot\mathbf{\hat{r}}_{ij}<0\\\
0,&\text{otherwise,}\end{cases}$ (32)

where $\rho_{ij}=(\rho_{i}+\rho_{j})/2$, and the signal velocity
$v_{\text{sig},i}$ is

$v_{\text{sig},i}=\alpha_{i}^{\text{AV}}c_{\text{s},i}+\beta^{\text{AV}}\left|\mathbf{v}_{ij}\cdot\mathbf{\hat{r}}_{ij}\right|,$
(33)

and $\alpha_{i}^{\text{AV}}\in[0,1]$ and $\beta^{\text{AV}}=2$. We use the
method of Morris & Monaghan (1997) to evolve $\alpha_{i}^{\text{AV}}$.
Finally, the contribution of the artificial viscosity to the gas internal
energy is

$\left(\derivative{u_{i}}{t}\right)^{\text{diss}}=-\frac{1}{\Omega_{i}}\sum_{j}\frac{m_{j}}{\rho_{ij}}v_{\text{sig},i}\frac{1}{2}(\mathbf{v}_{ij}\cdot\mathbf{\hat{r}}_{ij})^{2}F_{ij}(h_{i}).$
(34)

###  2.2 Implicit dust evolution algorithm

Equation 22 is in the explicit form for the evolution of the dust variable,
$s$. Here we present a method for solving the dust diffusion equation
implicitly using a backwards Euler method. This formalism uses Gauss-Seidel
iterations to solve the implicit dust evolution equation. The idea of using
this method to solve the dust evolution equation implicitly was inspired by
Whitehouse, Bate & Monaghan (2005), who used this method to solve the
diffusion equation for radiative transfer in the flux-limited diffusion
approximation. Using the implicit method detailed below allows us to use
timesteps that are not constrained by the stopping time of the dust, instead
using the hydrodynamical timestep and so typically speeding up calculations.

To advance a time dependent variable, $A$, from time $t=n$ to $t=n+1$ the
backwards Euler method can be stated as

$A_{i}^{(n+1)}=A_{i}^{(n)}+\mathrm{d}t\left(\derivative{A_{i}}{t}\right)^{(n+1)}.$
(35)

This gives us the implicit equation for $A$. For interactions between
particles $i$ and $j$, a Gauss-Seidel iterative method $A_{i}^{(n+1)}$ can be
solved for by rearranging the implicit equation to the form

$A_{i}^{(n+1)}=\frac{\displaystyle
A_{i}^{(n)}-\mathrm{d}t\sum_{j}\sigma_{ij}\left(A_{j}^{(n+1)}\right)}{\displaystyle
1-\mathrm{d}t\sum_{j}\sigma_{ij}},$ (36)

where $\sigma_{ij}$ contains quantities other than $A$. In this scheme
$A_{j}^{(n+1)}$ starts off as $A_{j}^{(n)}$ and is updated to $A_{j}^{(n+1)}$
as soon as this value becomes available. Eq. 36 is iterated over until some
specified convergence criterion is met. Taking Eq. 22, and putting it into
backwards Euler form we obtain

$\begin{split}s_{i}^{(n+1)}=s_{i}-&\mathrm{d}t\left\\{\frac{1+\left(s_{i}^{(n+1)}\right)^{2}}{2\rho_{i}}t_{\text{s}_{i}}\sum_{j}s_{j}^{(n+1)}L_{ij}+\right.\\\
&\left.\frac{\left(1+\left(s_{i}^{(n+1)}\right)^{2}\right)^{2}}{2\rho_{i}}\sum_{j}\frac{t_{\text{s}_{j}}}{1+\left(s_{j}^{(n+1)}\right)^{2}}s_{j}^{(n+1)}L_{ij}\right\\},\end{split}$
(37)

where
$L_{ij}\equiv\frac{m_{j}}{\rho_{j}}(P_{i}-P_{j})\frac{\overline{F}_{ij}}{r_{ij}}$.
The implicit equation for $s_{i}^{(n+1)}$ can be rearranged into a quartic
equation of the the form

$a_{4}\left(s_{i}^{(n+1)}\right)^{4}+a_{3}\left(s_{i}^{(n+1)}\right)^{3}+a_{2}\left(s_{i}^{(n+1)}\right)^{2}+a_{1}s_{i}^{(n+1)}+a_{0}=0,$
(38)

where

$\displaystyle a_{4}$
$\displaystyle=\frac{\mathrm{d}t}{2\rho_{i}}\sum_{j}\frac{t_{\text{s}_{j}}}{1+\left(s_{j}^{(n+1)}\right)^{2}}s_{j}^{(n+1)}L_{ij}^{(n+1)},$
$\displaystyle a_{3}$ $\displaystyle=0,$ $\displaystyle a_{2}$
$\displaystyle=\frac{\mathrm{d}t}{\rho_{i}}\sum_{j}\frac{t_{\text{s}_{j}}}{1+\left(s_{j}^{(n+1)}\right)^{2}}s_{j}^{(n+1)}L_{ij}^{(n+1)}+\mathrm{d}t\frac{t_{\text{s}_{i}}}{2\rho_{i}}\sum_{j}s_{j}^{(n+1)}L_{ij}^{(n+1)},$
$\displaystyle a_{1}$ $\displaystyle=1,$
$\displaystyle\begin{split}a_{0}&=\frac{\mathrm{d}t}{2\rho_{i}}\sum_{j}\frac{t_{\text{s}_{j}}}{1+\left(s_{j}^{(n+1)}\right)^{2}}s_{j}^{(n+1)}L_{ij}^{(n+1)}+\\\
&\quad\mathrm{d}t\frac{t_{\text{s}_{i}}}{2\rho_{i}}\sum_{j}s_{j}^{(n+1)}L_{ij}^{(n+1)}-s_{i}^{(n)}.\end{split}$

A solution to Eq. 38 can be found with any method of solving quartic
equations, e.g. numeric, however as $a_{\mathrm{3}}=0$ this is a depressed
quartic and we can solve the resolvent cubic (see Appendix A for the analytic
solution of the quartic). This analytic solution is preferred as it finds a
solution more quickly. Iterating the solution to Eq. 38 over all SPH particles
until convergence gives the values of the dust variable, $s$, at the next time
step. Again, we note that $s_{j}^{(n+1)}$ starts off as $s_{j}^{(n)}$ until we
have an updated value for $s_{j}$.

###  2.3 Convergence criterion

The value of $s_{i}^{(n+1)}$ is iterated over until convergence within an
accepted tolerance on the $m^{\text{th}}$ iteration, i.e. iterations continue
until

$\left|\frac{s_{i}^{n+1,m}-s_{i}^{n+1,m-1}}{s_{i}^{n+1,m}}\right|<\delta,$
(39)

where typically $\delta=10^{-3}$ and $s_{i}^{n+1,m}$ is the backward Euler
form of $s_{i}^{n+1}$ (equation 37) evaluated at iteration $m$. If the method
fails to converge, the timestep is split into two and the iterations begin
again. We continue in this manner until convergence is reached, or until the
original timestep has been split so many times that it is no longer
computationally viable.

In the limit where $\epsilon\geq 2.5\times 10^{-7}$ $(s>0.0005)$ the equation
is solved in this way. If the dust fraction drops below this value equation 38
is solved linearly (i.e. simply considering the last two terms on the left-
hand side) as the leading order terms of the quartic become vanishingly small.
We find no measurable loss of accuracy by solving the diffusion equation
linearly in this limit.

##  3 Numerical tests

We have implemented both the implicit and explicit forms of the dust-as-
mixture equations using the Ballabio et al. (2018) dust parameterisation into
the 3D SPH code sphNG. Originally sphNG was developed by Benz (1990) and Benz
et al. (1990) and subsequently modified significantly as described in Bate,
Bonnell & Price (1995), Whitehouse et al. (2005), Price & Monaghan (2007) and
other papers.

To test the dust-as-mixture dust fraction implementations, we used the
dustywave and dustshock tests as outlined in Laibe & Price (2011). In these
tests we adopt a stopping time defined as

$\displaystyle t_{\mathrm{s}}=\frac{\epsilon(1-\epsilon)\rho}{K},$ (40)

to be consistent with the literature and analytic solutions, where $K$ is the
coefficient of drag between the gas and the dust.

###  3.1 Dustywave

Figure 1: Results from the dustywave problem performed with our implicit
method with a varying drag coefficient $K=1,\ 10,\ 100,\ 1000$ (from top-left
to bottom-right) and a dust-to-gas ratio of 1:1. The analytic solution of the
gas velocity is plotted in red as a solid line and the analytic solution of
the dust velocity us plotted as a dashed line. The numerical solution of dust
and gas velocities are plotted at time $t=4.5$ as open and filled circles
respectively. As with an explicit method, the solution becomes inaccurate at
lower drag ($K\lesssim 40$) due to the breakdown of the terminal velocity
approximation.

The dustywave problem consists of a linear wave propagating through a dust-gas
mixture. The problem was originally described for dust-as-particles by Laibe &
Price (2011), then in dust-as-mixture by Laibe & Price (2014c). We set the
problem up in 3D, but only propagate the wave in the $x$ direction. The
problem is initialised with a mixture consisting of equal dust and gas
($\rho_{\mathrm{d}}=\rho_{\mathrm{g}}=1$, i.e. $\rho_{0}=2$, $\epsilon=0.5$)
with a sinusoidal perturbation in both density and velocity of the form
$v(x)=v_{0}\sin{(2\pi x)}$ and
$\rho(x)=\rho_{0}\left[1+\delta\rho_{0}\sin{(2\pi x)}\right]$. The amplitude
of the perturbation is set to $v_{0}=\delta\rho_{0}=1\times 10^{-4}$, with a
corresponding perturbation in the energy. An adiabatic equation of state is
used where $\gamma=5/3$ and the sound speed is set to $c_{\mathrm{s}}=1$. The
domain is $x,\ y,\ z\in[-0.52,0.52],\ [-0.069,0.069],\ [-0.073,0.073]$
respectively, in which 3744 particles are set up in a closed packed lattice,
giving a spacing of $\Delta p=0.02$, to ensure particle stability throughout
the calculation.

As noted by Price & Laibe (2015), there is an inconsistency between the
analytic and numerical solution. In the analytic solution
$\Delta\boldsymbol{v}$ is assumed to be zero at $t=0$, however in the terminal
velocity approximation $\Delta\boldsymbol{v}=t_{\mathrm{s}}\boldsymbol{f}$ is
non-zero. This inconsistency causes a small phase difference between the
analytic and numerical solution that gets larger as $t_{\mathrm{s}}$ gets
larger, i.e. less drag / larger grains, and specifically as the terminal
velocity approximation breaks down.

The solution is plotted after $\approx 4.5$ wave periods in Fig. 1 in which
both gas and dust velocities are plotted in filled and open circles
respectively. These values are reconstructed from the barycentric values for
velocity as follows

$\displaystyle\boldsymbol{v}_{\mathrm{g}}$
$\displaystyle=\boldsymbol{v}-\epsilon\Delta\boldsymbol{v},$ (41)
$\displaystyle\boldsymbol{v}_{\mathrm{d}}$
$\displaystyle=\boldsymbol{v}+(1-\epsilon)\Delta\boldsymbol{v}.$ (42)

The results from our implicit method are consistent with those of Price &
Laibe (2015) for the explicit method implemented using a direct second
derivative. For both explicit and implicit methods, the solutions are very
close to the analytic solutions in the regime where the terminal velocity
approximation is valid ($K\gtrsim 40$). In lower drag regimes, the terminal
velocity approximation gives a solution that is slightly out of phase or is
over-dampened in the lowest drag case (see the two upper panels of Fig. 1).

###  3.2 Dustyshock

We perform the dustyshock (Laibe & Price, 2012a) test here with our implicit
scheme at high drag. We use the standard Sod (1978) shock tube initial
conditions. In the gas for $x\leq 0$
$(\rho_{\mathrm{g}},v_{\mathrm{g}},P)=(1,0,1)$, and for $x>0$
$(\rho_{\mathrm{g}},v_{\mathrm{g}},P)=(0.125,0,0.125)$. We use an initial
dust-to-gas ratio of unity (i.e., $\epsilon=0.5$), and an adiabatic equation
of state where $\gamma=5/3$. The domain in 3D is $x,\ y,\ z\in[-0.5,0.5],\
[-0.012,0.012],\ [-0.017,0.017]$ containing 3509 particles, a particle spacing
of $\Delta p\approx 6.84\times 10^{-3}$. The numerical solution is shown in
Fig. 2 along with the analytic solution for the Sod shock tube in this
configuration. The numeric and analytic solutions are in good agreement, with
a slight discontinuity in the SPH solution at the contact discontinuity.

Figure 2: Results from the dustyshock problem with a high drag coefficient and
a high dust-to-gas ratio (1:1). The analytic solution is plotted at time
$t=0.2$ with a solid red line, the numerical solution is plotted with open
circles. From upper left going clockwise plotted are: velocity in $x$
direction, internal energy, gas pressure, and total density.

##  4 Dusty protostellar collapse

As the first application of this new algorithm for solving the dust diffusion
equations in dust-as-mixture SPH, we model the dynamics of small dust grains
during the protostellar collapse of a molecular cloud core. We compare the
results to those obtained using the semi-implicit dust-as-particles SPH method
(Lorén-Aguilar & Bate, 2014, 2015) to model the dynamics of the dust grains.
The dust-as-particles calculations are almost identical to those presented by
Bate & Loren-Aguilar (2017), except that the initial conditions are modified
to set the dust particles to their terminal velocity, since by construction
the dust-as-mixture SPH method discussed in this paper assumes the terminal
velocity approximation (Bate & Loren-Aguilar, 2017, took the dust particles to
be initially at rest relative to the gas). We not that this implementation of
dust-as-particles does not use the reconstruction of gas and dust velocities
of Price & Laibe (2020), so as a result linear waves present in the
calculation will tend to be overdamped when the dust particle separation
becomes too large to resolve the drag length-scale.

In addition to comparing the results from the dust-as-mixture and dust-as-
particles calculations, we perform the dust-as-mixture calculations using both
explicit and implicit methods for solving the dust diffusion equation to
compare their computational performance.

###  4.1 Initial conditions

Apart from the initial dust velocity being set to its terminal velocity, we
replicate the calculations presented in Bate & Loren-Aguilar (2017). We set up
unstable Bonnor-Ebert spheres of mass $5~{}\mathrm{M}_{\odot}$, radius $0.1$
pc, and an inner-to-outer density ratio of 20. We use spherical, reflective
boundary conditions. We use initially uniformly rotating clouds, with a
rotation rate of $1.012\times 10^{-13}~{}\mathrm{rad~{}s}^{-1}$, and a non-
rotating cloud. The magnitude of the ratio of rotational to gravitational
potential energy for the rotating cloud is $\beta=0.08$.

We use the same equation of state as Bate & Loren-Aguilar (2017) to model the
gas during different phases of collapse, as we do not use radiative transfer.
We use a barotropic equation of state for the gas pressure given by

$P=\left\\{\begin{array}[]{l
l}c_{\text{s,0}}^{2}\rho;&\rho<\rho_{\text{c}},\\\
c_{\text{s,0}}^{2}\rho_{\text{c}}\left(\rho/\rho_{\text{c}}\right)^{7/5};&\rho_{\text{c}}\leq\rho,\end{array}\right.$
(43)

where $c_{\text{s,0}}=1.87\times 10^{4}$ cm s-1 is the initial isothermal
sound speed of gas with a temperature of 10 K (the mean molecular weight is
$\mu=2.38$), and $\rho_{\text{c}}=10^{-13}$ g cm-3. As in Bate & Loren-
Aguilar, our calculations only follow the collapse until shortly after the
density exceeds $\rho_{\rm c}$, so the gas almost always isothermal at 10 K,
except in the centre at the end of the calculations when the temperature
begins to rise.

In the dust-as-particles calculations we use 1 million SPH particles to model
the gas, and 300,000 to model the dust. Comparatively in the dust-as-mixture
calculations we use 1 million SPH particles that model both gas and dust as
one fluid in a dust-as-mixture model. A uniform initial dust-to-gas ratio of
1:100 is used for all models. Using identical numbers of SPH particles to
model the gas components in both sets of calculations means the gas
resolutions are the same and the only difference between the calculations is
due to the dust implementation.

###  4.2 Results

Figure 3: Time evolution of the central density of gas and dust during the
collapse of non-rotating molecular cloud cores performed using in dust-as-
mixture (left) and dust-as-particles (right) methods. There is a monatonic
increase in the maximum density of dust with time. The dust-as-mixture and
dust-as-particles solutions are similar, except for the largest dust grains.
In the dust-as-particles $100\mu\mathrm{m}$ calculation, the dust in the outer
parts of the cloud is collapses more rapidly than in the inner regions,
resulting in a ‘pile up’ of dust. For this case, the terminal velocity
approximation of the dust-as-mixture method is only comparable until $t\approx
70,000$ yrs, so we only plot the dust-as-particles solution to this point. The
gas density is from the $0.1\mu\mathrm{m}$ calculations.

Figure 4: Time evolution of the central density of gas and dust in both
rotating calculations in dust-as-mixture (left) and dust-as-particles (right).
Due to the rotation of the cloud the formation of the hydrostatic core is
delayed. The main differences between the two methods is with the largest
grains; it is with this size of grain where the underlying assumptions of
dust-as-mixture begin to breakdown and a dust-as-particles solution is
desirable. The gas density is from the $0.1\mu\mathrm{m}$ calculations. Figure
5: Evolution of dust grain column density viewed in the plane of rotation for
calculations of molecular cloud core collapse using the dust-as-mixture
method. Plotted from left to right are calculations with dust grain sizes
$r_{\mathrm{s}}=1\mu\mathrm{m},\ 30\mu\mathrm{m},\ 100\mu\mathrm{m}$. The
colour bar is set such that it is the same in each panel of the plot. The
small dust grains are coupled to the gas and follow the gas well. The larger
grains exhibit different behaviour whereby they collapse more quickly
resulting in the central oblate region having a higher dust density. The
central obliquity is due to the overall rotation of the gas cloud.  Figure 6:
Evolution of dust grain column density in the plane of cloud rotation for
calculations of molecular cloud core collapse using the dust-as-particles
method. Shown from left to right are calculations with dust grain sizes
$r_{\mathrm{s}}=1\mu\mathrm{m},\ 30\mu\mathrm{m},\ 100\mu\mathrm{m}$. The
colour bar is set such that it is the same in each panel of the plot. In the
$100\mu\mathrm{m}$ calculation the dust grains from the outer regions of the
cloud collapse faster than the rest of the grains causing a pile up of grains
that pass through each other in the centre of the cloud. These subsequently
settle into a dense flat disc, whereupon the calculation grinds to a halt.

Here we present the results of the protostellar collapse of a non-rotating
cloud and a rotating cloud with $\beta=0.08$ using both the dust-as-mixture
and -particles SPH methods.

In Fig. 3 we show the time evolution of the densities of gas and dust at the
centre of the collapsing non-rotating clouds. In the dust-as-mixture
calculations (left panel) there is a monotonic increase in the density of
dust, and the rate of increase is greater for larger dust grain sizes due to
the weaker coupling to the gas, although there is little difference for grains
of size $r_{\mathrm{s}}\lesssim 10\mu\text{m}$. The results are very similar
in the equivalent dust-as-particles calculations (right panels), except for
the largest grain size. For very small dust grain sizes there is good
agreement between the dust-as-mixture and dust-as-particles results that only
starts to break down very late in the collapse at very high densities as the
first hydrostatic core (Larson, 1969) forms; in the the dust-as-particles
calculations (unphysical) dust particle clumping occurs late in the collapse.
For intermediate dust grain sizes ($r_{\mathrm{s}}\approx 10-30~{}\mu$m) the
central dust density tends to be slightly higher for the dust-as-particles
calculations than the dust-as-mixture calculations in the latter third of the
collapse, again due to slight dust clumping when using the dust-as-particles
method.

In the dust-as-particles calculation with dust grains of size
$r_{\mathrm{s}}=100\mu\mathrm{m}$, the initial terminal infall velocities of
the dust in the outer regions of the cloud are much higher than in the centre.
Therefore, the dust in the outer regions collapses much faster than the dust
in dense inner regions causing a ‘pile-up’ as the dust collapses. This is an
extreme case for which the terminal velocity approximation assumed by the
dust-as-mixture model quickly breaks down. Therefore, for this case, we
present the dust-as-particles results in Fig. 3 only until the dust from the
outer regions of the cloud catches up with that from the inner part of the
cloud. The results between the dust-as-mixture and dust-as-particles for
$100\mu\mathrm{m}$ grains are similar until this point ($t\approx 70,000$
yrs).

In Fig. 4 we show the time evolution of the gas and dust densities at the
centre of the rotating clouds. In the calculations with dust grains
$r_{\mathrm{s}}<100\mu\mathrm{m}$ there is little difference in the central
densities of dust and gas between the two methods. Again, the dust-as-
particles calculations are prone to some dust clumping. Also again, the
largest grains ($100\mu\text{m}$) are less well coupled to the gas and so
collapse much more quickly than the small grains. In the rotating case the
collapse happens slower than in the non-rotating case, as expected due to
centrifugal forces. But again, the dynamics of the large dust grains is better
captured using the dust-as-particles method as in this case the terminal
velocity approximation is not appropriate.

In Fig. 5 we show snapshots of the evolution of the column density of dust in
the rotating protostellar collapse calculations containing dust with grain
sizes of $r_{\mathrm{s}}=1\mu\mathrm{m},\ 30\mu\mathrm{m},\ 100\mu\mathrm{m}$.
The region shown is the centre of the cloud in the plane of rotation. The
small dust grains are tightly coupled to the gas (the morphology of the gas
column density looks almost identical to the 1$\mu$m dust morphology), whereas
the larger grains are not as well coupled and collapse more quickly. This
leads to an enhancement in dust-to-gas ratio when larger dust grains are
present. The $100\mu\mathrm{m}$ dust grains collapse fast enough to form a
disc-like structure well before the stellar core forms, a consequence of the
rotation rate of the cloud. A qualitatively similar result was found by Bate &
Loren-Aguilar (2017), shown in Fig. 3 of their work. We show the equivalent
plot of Fig. 5 for the dust-as-particles calculations in Fig. 6. The dust
density for grain sizes 1 and 30$\mu\mathrm{m}$ are similar in both the one-
and dust-as-particles calculations, with the exception that in the dust-as-
particles calculation the density at the centre of the cloud is slightly
denser than in the dust-as-mixture at very late times. The 100$\mu\mathrm{m}$
dust grain calculation using dust-as-particles looks vastly different. In this
calculation the dust grains in the outer parts of the cloud are very decoupled
from the gas and have such a large initial terminal velocities that they are
able to collapse to the centre of the cloud very quickly. Moreover, because
the gas drag on them is so weak, even when they get to the denser inner
regions the drag is insufficient to slow them to their local terminal velocity
so they continue to collapse more quickly than in the dust-as-mixture
calculation. This is a case that can not be modelled accurately using the
dust-as-mixture terminal velocity approximation. We note that this is an
extreme case – such large dust grains are not thought to be abundant in
molecular cloud cores. Interstellar dust grains are typically thought to be
sub-micron in size.

###  4.3 Comparison with explicit timestepping

In this section we compare the computational performance of the explicit
timestepping algorithm of the dust evolution equation as described in Ballabio
et al. (2018) with the implicit algorithm described in this work. In their
formulation Ballabio et al. (2018) impose a criterion upon the stopping time
of the dust to prevent numerical artefacts in regions of low dust-to-gas
ratio, namely

$\tilde{t}_{\mathrm{s}}=\min(t_{\mathrm{s}},h/c_{\mathrm{s}}).$ (44)

The reasoning for adopting this stopping time is to prevent infeasibly slow
calculations by limiting the mass flux of particles in regions where there are
steep dust density gradients, i.e. at the edge of discs. However this
artificially alters the drag and thus the effective size of the dust grains
and should be used with caution.

As an example of the problems that altering the stopping time can give, we
repeat some of the dust-as-mixture non-rotating protostellar collapse
calculations from Section 4 using explicit time stepping with the stopping
time limiter turned on and off. We compare the results and performance with
our implicit algorithm without a stopping time limiter.

In Fig. 7 we show the distribution of dust as a function of radius in the
collapsing cloud when the maximum gas density reaches $\rho=10^{-10}\
\mathrm{g~{}cm}^{-3}$ using both the explicit and implicit methods. The
explicit results have the stopping time limiter applied, but no stopping time
limiter is used for the implicit calculations. We give results for dust grain
sizes $r_{\mathrm{s}}=0.1,\ 100\mu\mathrm{m}$. In the explicit, stopping time
limited calculations the dust density distribution for both grain sizes is
very similar throughout the cloud, and identical in the outer regions where
the gas density is low and the stopping time is the longest (i.e. where the
stopping time limiter is most likely to be active). The implicit result for
the $0.1\mu$m grain size is very similar to the explicit result. However, in
the implicit calculation with large grains the dust density is very different.
In the outer parts of the cloud, the dust has fallen inwards leaving no dust
in the outer parts of the cloud and producing much higher dust densities in
the inner regions. Near the centre of the clouds the dust density
distributions for the two explicit calculations with the stopping time limiter
differ slightly, with the larger grains falling faster leading to a somewhat
higher density than the small grains, due to the stopping time limited turning
off at the higher gas densities. However the density of the larger grains here
is still substantially underestimated compared with the implicit calculation.
The large grains in the explicit calculation behave as if they were much
smaller due to the timestep limiter. This is why extreme caution is needed if
a stopping time limiter is used – a plausible result is obtained for the
$r_{\mathrm{s}}=100\mu\mathrm{m}$ explicit calculation with the stopping time
limiter, but it is a physically incorrect result.

Figure 7: Dust density as a function of radius for dust grain sizes $0.1$ and
$100\mu\mathrm{m}$ in non-rotating calculations of protostellar collapse when
the peak gas density reaches of $10^{-10}\text{g cm}^{-3}$, just after the
first hydrostatic core begins to form. The ‘expl’ lines denote the dust
density from explicit dust-as-mixture calculations with the stopping time
limited applied, and ‘impl’ lines denote the equivalent calculations from
implicit dust-as-mixture calculations (without a stopping time limiter). Due
to the action of the stopping time limiter the dust density profiles for the
two explicit calculations with different grain sizes are very similar, whereas
with the implicit method the large grains have evacuated the outer parts of
the cloud and migrated much more strongly inward.

To compare the computational performance between the explicit algorithm and
the implicit algorithm, we also ran an explicit calculation without any
stopping time limiter for the $r_{\mathrm{s}}=0.1\mu\mathrm{m}$ dust grains (a
typical size for interstellar dust grains). This increased the calculation
time so significantly that it was unrealistic to allow the gas density to
reach $\rho=10^{-10}\mathrm{g~{}cm}^{-3}$. We estimate that the implicit
algorithm is at least two orders of magnitude (i.e., several hundred times)
faster than the explicit algorithm at calculating this protostellar collapse
case, because the implicit calculation takes standard hydrodynamical timesteps
and does not require the use of a dust-as-mixture dust timestep. When
comparing the dust density distributions of these two calculations they are
almost identical; our algorithm yields the same result as the explicit
algorithm without limiting the stopping time, but does so much more quickly.
When using the explicit method with the stopping time limiter for the
$0.1\mu$m grain size, the explicit calculation was only about 5 percent faster
than the implicit method.

##  5 Protoplanetary discs

In this section we present the results from protoplanetary disc simulations
similar to those of Ballabio et al. (2018). The calculations include an
embedded planet to study the diffusivity effects of gradients in dust fraction
introduced by planet-disc interactions. The idea is that as the planet orbits,
gas and dust are flung out into the outer regions of the disc where there are
steep dust fraction gradients. It is in these environments that Ballabio et
al. found that numerical artifacts in the dust can arise. This is most notable
when using earlier formulations that evolved the dust fraction or the dust-to-
gas ratio.

For our planet-in-disc test calculations, the initial radial extend of the gas
is $r\in[25,200]$ au and the dust initially extends to 90% of the gas radius.
A locally-isothermal equation of state is used, with a radial temperature
proportional to $r^{-1/2}$ so that ratio of the disc scale height to radius is
a constant, $H/r=0.1$. The central object is modelled using a gravitational
potential of mass $1~{}\mathrm{M}_{\odot}$. The surface density of the disc
has an initial profile proportional to $r^{-1/2}$. The mass of the disc is
$0.0348~{}\mathrm{M}_{\odot}$ and the dust-to-gas ratio is set at 0.01. We
used 1 million SPH particles to model the disc. The disc contains a 6 Jupiter-
mass (MJ) planet at a radius of $140$ au. The planet is modelled using a sink
particle (Bate et al., 1995), with an accretion radius of 10 au. The sink
particle only interacts with the gas via gravity and accretion, and the
trajectory of the sink particle is integrated with the same Runge-Kutta-
Fehlberg integrator that is used for the SPH particles, but with a much lower
tolerance.

Figure 8: The time evolution of the column density of dust in a protoplanetary
disc containing a 6 MJ planet. Upper panels: Results using the implicit dust-
as-mixture dust algorithm. Lower panels: Results using the explicit dust-as-
mixture dust algorithm with the stopping time limiter. Both calculations
contain $100\mu\mathrm{m}$ dust grains and use the same initial conditions.
Similar results are obtained, but the implicit method gives slightly ‘sharper’
edges to the dust in the gap produced by the planet and at the outer edge of
the disc.

In Fig. 8 we show the time evolution of the dust column density within the
protoplanetary disc containing dust grains of size $100\mu\mathrm{m}$. The
upper panels show the results from a calculation using the implicit dust-as-
mixture algorithm to evolve the dust, the low panels show an identical
calculation but using the explicit dust-as-mixture algorithm with the stopping
time limiter.

The overall dust density distributions are almost identical. However, there
are some differences in regions of high dust density gradient at the dust gap
created by the planet and at the outer edge of the dust disc. In the explicit
dust-as-mixture calculation the dust density is more smoothed out, whereas in
the implicit dust-as-mixture calculation the edges of the dust disc are more
sharply defined. This is most easily spotted in the gap in the dust; the dust
does not diffuse into the gap as much in the implicit calculation.

We tested the mass conservation of the implicit dust-as-mixture algorithm and
find it performs similarly to the explicit stopping time limited algorithm at
conserving mass. We performed a similar calculation to the one described above
except that we used only 200,000 SPH particles and we modelled the planet and
disc in the same way as Ayliffe & Bate (2009) where gas and dust falls on to
and is bound to a planetary surface and remains in the calculation (rather
than being accreted by a sink particle). Since all of the mass remains in the
calculation (as opposed to being accreted by the sink particle) this allows us
to measure more precisely the mass conservation of the respective dust
schemes. In the upper panel of Fig. 9 we show that we recover the result
reported by Ballabio et al. (2018) in which the parameterisation of dust
fraction given by Price & Laibe (2015), $\epsilon=s^{2}/\rho$, does not
conserve mass (the dust mass increases with time). In contrast, our implicit
dust algorithm maintains the dust mass as well as the explicit algorithm with
the Ballabio et al. (2018) dust parameterisation and stopping time limiter. In
the lower panel of Fig. 9 we focus on just the implicit and explicit
calculations to compare how well they converse mass against each other. Our
implicit calculation slowly loses a small fraction of mass ($\approx$ 0.02 %)
initially and then approximately maintains this mass throughout the remainder
of the calculation. In the case of the explicit calculation, the initial dust
mass very slightly decreases through most of the calculation, but then loses a
small fraction of mass towards the end ($\approx$ 0.06%). Overall we find that
the dust mass conservation of the two different algorithms is similarly good.

Figure 9: Dust mass evolution in a disc with an embedded planet. Upper: Dust
mass evolution for the implicit and explicit calculations using the Ballabio
et al. (2018) dust variable, and an explicit calculation using an alternate
dust variable. Lower: A zoom in of the results for our implicit dust evolution
algorithm and the explicit dust evolution algorithm, both using the Ballabio
et al. parameterisation. Both of these calculations conserve dust mass to
better than 0.1 percent.

##  6 Negativity of the dust parameter

Whilst the parameterisation of the dust fraction as given by in equation 1
ensures that it remains positive and less than one, there is no such
constraint on the dust evolution variable, $s$, itself. During calculations,
as the parameter is evolved, it can become negative. The dust fraction itself
remains limited to between zero and unity, but a negative dust variable can
still produce unwanted, nonphysical effects in the dust fraction. These
effects typically manifest as dust seemingly ‘leaking out’ into regions of low
dust fraction, at the edges of discs for example. This was pointed out by
Price & Laibe (2015) for earlier versions of the dust parameterisation. They
found that the most effective way to counter this issue was to set the dust
fraction of the problem particles to zero.

In what follows, we show that the problem of the dust variable going negative
can still occur with the newer dust parameterisation of Ballabio et al.
(2018). Here we demonstrate the effects of a negative dust variable in both an
idealised test case and in a more physically motivated fully 3D simulation. We
then demonstrate that another advantage of the implicit dust method is that it
avoids this problem entirely.

###  6.1 Dust settling test

For the idealised test case, we use the dust settling test of Price & Laibe
(2015) which considers vertical dust settling in a Cartesian box with an
acceleration in the ‘vertical’ (z) direction that mimics a small piece of a
protoplanetary disc (i.e., the isothermal gas has a Gaussian density
distribution vertically centred on $z=0$). A full description of the set up
the test is given by Price & Laibe (2015) and will not be repeated here. In
our case, we start the test with a dust-to-gas ratio of 1:100 and use a dust
grain size of 100$\mu$m. The orbital period is 353 years and we show the state
of the dust layer after 20 orbits.

As the disc evolves, the dust falls out of the low-density surface regions of
the disc to create a dust layer with steep dust gradients at its surface. We
show our results in Fig. 10 using both the explicit and implicit methods. The
problem of the dust variable becoming negative occurs at this sharp gradient
at the ‘surface’ of the dust layer. This is due to the way exact mass
conservation of the dust-as-mixture equations is ensured by pairwise SPH
particle interactions. The equations ‘take’ dust from one SPH particle and
‘give’ it to another. However, if the dust variable of the donor particle is
small this can lead to the dust variable of the donor particle becoming
negative. Due to the definition of the dust fraction, a negative dust variable
still gives a positive dust fraction. However, this means that in regions
where no dust is expected, dust diffuses out into these regions in an
unphysical manner (lower panel of Fig. 10). This can be avoided in the
explicit method by explicitly forcing the dust variable to remain positive,
but this comes at the expense of sacrificing dust mass conservation.

Using our implicit method, in which the solution to the dust diffusion
equation is iterated using the Gauss-Seidel method, the dust variable remains
positive and the apparent diffusion of the dust out of the dust layer is
avoided (upper panel of Fig. 10).

Figure 10: The dust settling test with 100$\mu$m dust grains after 20 orbital
periods. Upper panel: dust fraction as a function of vertical height in the
disc using the implicit method to evolve the dust. Lower panel: the same
except using the explicit method. Using the explicit method, dust ‘leaks’ out
of the dust layer into the low-density atmosphere of the disc ($|z|>6$ au in
the above panels). Using the implicit method gives the same dust profile near
the midplane of the disc, but outside the dust layer where the dust fraction
should be zero, it is zero as it should be.

###  6.2 Protoplanetary disc

Here we demonstrate the effect of the dust variable going negative in regions
with a large dust fraction gradient also occurs in more complex simulations.
From the same protoplanetary disc calculations that were discussed in Section
5 we take a cross section through both discs. This is shown in Fig. 11, the
upper panel shows the dust density cross section of the disc from the implicit
calculation and the lower is the same for the explicit algorithm. Whilst the
bulk of the dust across both calculations is settling towards the midplane of
the disc, in the calculation that used the explicit dust algorithm there is an
apparent excess of dust above and below the the bulk of the dust layer. In
addition to this enhancement above and below the midplane, there is also
excess dust in the gap carved by the planet. This is due to the dust variable
going negative in these regions. This is the same effect as that demonstrated
by the simple dust settling test but in a full 3D global simulation. Again,
the implicit method for evolving the dust diffusion avoids this problem.

Figure 11: Cross section of the dust density in the protoplanetary disc with
an embedded planet. Upper panel: Calculation using the implicit dust
algorithm. Lower panel: Calculation using the explicit dust algorithm. In the
explicit calculation there is an enhancement in dust density above (and below)
the bulk of where the dust is expected to be. This is an artefact of the dust
variable becoming negative. Using the implicit method avoids this problem.

##  7 Conclusions

In this paper we have described an implicit algorithm to solve the dust-as-
mixture dust diffusion equation of Price & Laibe (2015). We have benchmarked
this algorithm against standard tests in the literature. The algorithm has
also been tested in the settings of the collapse of a dense molecular cloud
core and in a protoplanetary disc containing a planet, and compared with
results and methods in the literature. The main features of the algorithm we
have described are summarised as follows

  1. 1.

The timestep criterion required by explicit dust diffusion can be avoided
meaning that the dust evolution can be evolved on hydrodynamical timesteps
which, depending on the dust size, can be much larger.

  2. 2.

The time to solution compared to a fully explicit solution can be greatly
reduced in the regimes of intermediate to large dust grains (by orders of
magnitude). Even for very small dust grains, the implicit method is not much
slower.

  3. 3.

Due to the above there is no need to artificially limit the stopping time of
dust particles, which has the adverse effect of essentially changing the dust
grain size.

  4. 4.

The iterative solution of the dust diffusion equation in the implicit method
naturally avoids the dust variable becoming negative. This avoids dust
apparently ‘leaking’ through strong dust fraction gradients.

Overall we find this implicit algorithm to be fast computationally, to solve
the dust diffusion equations accurately, and to avoid several problems that
manifest in explicit calculations. However, as with any dust-as-mixture
method, the method still fails to give accurate solutions in the limit of
large (weakly coupled) grains for which the stopping time is long. Using a
dust-as-particles method is more appropriate in such circumstances.

## Acknowledgements

The authors thank the anonymous referee for their useful comments. DE thanks
Pablo Loren-Aguilar for useful discussions, and Thomas J. R. Bending for
guidance regarding the sphNG code. DE is funded by a Science and Technology
Facilities Council (STFC) studentship (ST/V506679/1). The calculations
discussed in this paper were performed on the University of Exeter
Supercomputer, ISCA, and on the DiRAC Data Intensive system, operated by the
University of Leicester IT Services, which forms part of the STFC DiRAC HPC
Facility (www.dirac.ac.uk). The latter equipment is funded by BIS National
E-Infrastructure capital grant ST/K000373/1 and STFC DiRAC Operations grant
ST/K0003259/1. DiRAC is part of the National e-Infrastructure. Several figures
were made using splash (Price, 2007) and matplotlib (Hunter, 2007).

## Data Availability

The data and scripts used to produce the figures in this article can be found
at https://doi.org/10.5281/zenodo.10794912. A description of the files and
scripts are also provided.

## References

  * Abramowitz & Stegun (1972) Abramowitz M., Stegun I. A., 1972, Handbook of mathematical functions. US Government printing office
  * Ayliffe & Bate (2009) Ayliffe B. A., Bate M. R., 2009, Monthly Notices of the Royal Astronomical Society, 393, 49
  * Bai & Stone (2010) Bai X.-N., Stone J. M., 2010, The Astrophysical Journal Supplement Series, 190, 297
  * Ballabio et al. (2018) Ballabio G., Dipierro G., Veronesi B., Lodato G., Hutchison M., Laibe G., Price D. J., 2018, MNRAS, 477, 2766
  * Barrière-Fouchet et al. (2005) Barrière-Fouchet L., Gonzalez J. F., Murray J. R., Humble R. J., Maddison S. T., 2005, Astronomy and Astrophysics, 443, 185
  * Bate & Loren-Aguilar (2017) Bate M. R., Loren-Aguilar P., 2017, Monthly Notices of the Royal Astronomical Society, 465, 1089
  * Bate et al. (1995) Bate M. R., Bonnell I. A., Price N. M., 1995, MNRAS, 277, 362
  * Benz (1990) Benz W., 1990, in Buchler J. R., ed., Numerical Modelling of Nonlinear Stellar Pulsations Problems and Prospects. p. 269
  * Benz et al. (1990) Benz W., Bowers R. L., Cameron A. G. W., Press W. H. ., 1990, ApJ, 348, 647
  * Benítez-Llambay et al. (2019) Benítez-Llambay P., Krapp L., Pessah M. E., 2019, The Astrophysical Journal Supplement Series, 241, 25
  * Hunter (2007) Hunter J. D., 2007, Computing in Science & Engineering, 9, 90
  * Hutchison et al. (2016) Hutchison M. A., Price D. J., Laibe G., Maddison S. T., 2016, Monthly Notices of the Royal Astronomical Society, 461, 742
  * Hutchison et al. (2018) Hutchison M., Price D. J., Laibe G., 2018, Monthly Notices of the Royal Astronomical Society, 476, 2186
  * Hutchison et al. (2022) Hutchison M. A., Bodénan J.-D., Mayer L., Schönbächler M., 2022, Monthly Notices of the Royal Astronomical Society, 512, 5874
  * Laibe & Price (2011) Laibe G., Price D. J., 2011, Monthly Notices of the Royal Astronomical Society, 418, 1491
  * Laibe & Price (2012a) Laibe G., Price D. J., 2012a, Monthly Notices of the Royal Astronomical Society, 420, 2345
  * Laibe & Price (2012b) Laibe G., Price D. J., 2012b, Monthly Notices of the Royal Astronomical Society, 420, 2365
  * Laibe & Price (2014a) Laibe G., Price D. J., 2014a, Monthly Notices of the Royal Astronomical Society, 440, 2136
  * Laibe & Price (2014b) Laibe G., Price D. J., 2014b, Monthly Notices of the Royal Astronomical Society, 440, 2147
  * Laibe & Price (2014c) Laibe G., Price D. J., 2014c, Monthly Notices of the Royal Astronomical Society, 444, 1940
  * Larson (1969) Larson R. B., 1969, MNRAS, 145, 271
  * Lebreuilly et al. (2019) Lebreuilly U., Commerçon B., Laibe G., 2019, Astronomy and Astrophysics, 626, A96
  * Lorén-Aguilar & Bate (2014) Lorén-Aguilar P., Bate M. R., 2014, Monthly Notices of the Royal Astronomical Society, 443, 927
  * Lorén-Aguilar & Bate (2015) Lorén-Aguilar P., Bate M. R., 2015, Monthly Notices of the Royal Astronomical Society, 454, 4114
  * Maddison et al. (2003) Maddison S. T., Humble R. J., Murray J. R., 2003, in Deming D., Seager S., eds, Astronomical Society of the Pacific Conference Series Vol. 294, Scientific Frontiers in Research on Extrasolar Planets. pp 307–310 (arXiv:astro-ph/0209153), doi:10.48550/arXiv.astro-ph/0209153
  * Monaghan (1997) Monaghan J. J., 1997, Journal of Computational Physics, 138, 801
  * Monaghan & Kocharyan (1995) Monaghan J. J., Kocharyan A., 1995, Computer Physics Communications, 87, 225
  * Morris & Monaghan (1997) Morris J. P., Monaghan J. J., 1997, Journal of Computational Physics, 136, 41
  * Price (2007) Price D. J., 2007, Publications of the Astronomical Society of Australia, 24, 159
  * Price & Laibe (2015) Price D. J., Laibe G., 2015, MNRAS, 451, 813
  * Price & Laibe (2020) Price D. J., Laibe G., 2020, Monthly Notices of the Royal Astronomical Society, 495, 3929
  * Price & Monaghan (2007) Price D. J., Monaghan J. J., 2007, MNRAS, 374, 1347
  * Sod (1978) Sod G. A., 1978, Journal of Computational Physics, 27, 1
  * Stoyanovskaya et al. (2018) Stoyanovskaya O. P., Glushko T. A., Snytnikov N. V., Snytnikov V. N., 2018, Astronomy and Computing, 25, 25
  * Stoyanovskaya et al. (2020) Stoyanovskaya O. P., Okladnikov F. A., Vorobyov E. I., Pavlyuchenkov Y. N., Akimkin V. V., 2020, Astronomy Reports, 64, 107
  * Whitehouse et al. (2005) Whitehouse S. C., Bate M. R., Monaghan J. J., 2005, MNRAS, 364, 1367
  * Youdin & Goodman (2005) Youdin A. N., Goodman J., 2005, ApJ, 620, 459

##  Appendix A Solving the quartic analytically

The solutions to the quartic equation

$P(x)=x^{4}+a_{3}x^{3}+a_{2}x^{2}+a_{1}x+a_{0}=0,$ (45)

can be found via solving the resolvent cubic if $P(x)$ is depressed, i.e.
$a_{3}=0$. Firstly, the real root, $y_{1}$ is found for the resolvent cubic

$y^{3}-a_{2}y^{2}+-4a_{0}y+4a_{2}a_{0}-a_{1}^{2}=0.$ (46)

Then the four solutions to the quartic $P(x)$ are given by the four solutions
of (Abramowitz & Stegun, 1972)

$z^{2}\mp\left(y_{1}-a_{2}\right)^{1/2}z+\left(\frac{y_{1}}{2}\mp\left(\frac{y_{1}}{2}\right)^{2}-a_{0}\right)=0.$
(47)


# Some extended hypergeometric matrix functions and their fractional calculus

Ashish Verma 111Department of Mathematics, Prof. Rajendra Singh (Rajju
Bhaiya), Institute of Physical Sciences for Study and Research V.B.S.
Purvanchal University, Jaunpur (U.P.)- 222003, India.
Email<EMAIL_ADDRESS>(Corresponding author), Ravi Dwivedi 222Department
of Basic Science and Humanities, PSIT College of Engineering, Kanpur 209305,
India.
Email<EMAIL_ADDRESS>and Vivek Sahai 333Department of Mathematics and
Astronomy, Lucknow University, Lucknow 226007, India.
Email<EMAIL_ADDRESS>
###### Abstract

In this paper, we study some extended hypergeometric functions from matrix
point of view. We have given the integral representations of these matrix
functions. Finally, we obtain some generating function relations using
fractional derivative operators.

Keywords: Matrix functional calculus, Appell functions, Lauricella functions,
Fractional derivative operators

AMS Subject Classification: 15A15; 33C65.

##  1 Introduction

The Gauss hypergeometric function and Kummer hypergeometric function are
amongst the most appeared special functions in mathematics as well as in
physics. Due to a wide range of applications, these hypergeometric functions
have been studied from different aspects _viz._ from matrix point of view [10,
12, 14]; extended Gauss and Kummer hypergeometric functions [4] and finite
field analogue [9] to name a few. Recently, the matrix version of extended
beta function, extended Gauss hypergeometric function and Kummer
hypergeometric function have been studied by Abdalla and Bakhet in [1, 2]. The
regions of convergence, integral representations and differential formulas
satisfied by these matrix functions are determined. Verma discussed recursion
formula, infinite summation formula for the Srivastava’s triple hypergeometric
matrix functions $H_{\mathcal{A}}$, $H_{\mathcal{B}}$ and $H_{\mathcal{C}}$
[18]. In [15, 16, 17, 19], authors introduced the incomplete first, second and
fourth Appell hypergeometric matrix functions and incomplete Srivastava’s
triple hypergeometric matrix functions and studied some basic properties;
matrix differential equation, integral formula, recursion formula, recurrence
relation and differentiation formula of these functions. In this paper, we
study the matrix analogues of some extended hypergeometric functions of two
and three variables. We obtain the integral representations for these matrix
functions. Finally, we use the extended Riemann-Liouville derivative operator
and obtain certain generating function relations in terms of these matrix
functions. The section-wise treatment is as follows:

In Section 2, we give the basic definitions related to special matrix
functions that are needed in the sequel. In Section 3, we define the extended
Appell matrix functions and extended Lauricella matrix function. We also
discuss the regions of convergence as well as the integral representations
satisfied by these matrix functions. In Section 4, we define the matrix
analogue of extended Riemann-Liouville derivative operators and transform
several matrix functions under these derivative operators. Finally, in Section
5, we obtain certain generating function relations which turn out in terms of
extended Appell matrix functions and extended Lauricella matrix function.

##  2 Preliminaries

Throughout the paper, $\mathbb{C}^{r\times r}$ is the vector space of
$r$-square matrices with complex entries. For a matrix
$A\in\mathbb{C}^{r\times r}$ the spectrum, denoted by $\sigma(A)$, is the set
of eigenvalues of $A$. If $\Re(z)$ denotes the real part of a complex number
$z$, then a matrix $A$ in $\mathbb{C}^{r\times r}$ is said to be positive
stable if $\Re(\lambda)>0$ for all $\lambda\in\sigma(A)$.

If $A$ is a positive stable matrix in $\mathbb{C}^{r\times r}$, then
$\Gamma(A)$ can be expressed as [10]

$\Gamma(A)=\int_{0}^{\infty}e^{-t}\,t^{A-I}\,dt.$ (2.1)

Furthermore, if $A+kI$ is invertible for all integers $k\geq 0$, then the
reciprocal gamma matrix function is defined as [10]

$\Gamma^{-1}(A)=A(A+I)\dots(A+(n-1)I)\Gamma^{-1}(A+nI),\ n\geq 1.$ (2.2)

If $M\in\mathbb{C}^{r\times r}$ is a positive stable matrix and $n\geq 1$ is
an integer, then the gamma matrix function can also be defined in the form of
a limit as [10]

$\Gamma(M)=\lim_{n\to\infty}(n-1)!\,(M)_{n}^{-1}\,n^{M}.$ (2.3)

By application of the matrix functional calculus, the Pochhammer symbol for
$A\in\mathbb{C}^{r\times r}$ is given by [11]

$(A)_{n}=\begin{cases}I,&\text{if $n=0$},\\\ A(A+I)\dots(A+(n-1)I),&\text{if
$n\geq 1$}.\end{cases}$ (2.4)

This gives

$(A)_{n}=\Gamma^{-1}(A)\ \Gamma(A+nI),\qquad n\geq 1.$ (2.5)

If $A$ and $B$ are positive stable matrices in $\mathbb{C}^{r\times r}$, then,
for $AB=BA$, the beta matrix function is defined as [10]

$\displaystyle\mathfrak{B}(A,B)=\Gamma(A)\Gamma(B)\Gamma^{-1}(A+B)$
$\displaystyle=\int_{0}^{1}t^{A-I}(1-t)^{B-I}dt$ (2.6)
$\displaystyle=\int_{0}^{\infty}u^{A-I}(1+u)^{-(A+B)}du.$ (2.7)

Let $A$, $B$ and $\mathbb{X}$ be positive stable and commuting matrices in
$\mathbb{C}^{r\times r}$ such that $A+kI$, $B+kI$ and $\mathbb{X}+kI$ are
invertible for all integer $k\geq 0$. Then the extended beta matrix function
$\mathfrak{B}(A,B;\mathbb{X})$ is defined by [1]

$\displaystyle\mathfrak{B}(A,B;\mathbb{X})=\int_{0}^{1}t^{A-I}(1-t)^{B-I}\exp\Big{(}\frac{-\mathbb{X}}{t(1-t)}\Big{)}dt.$
(2.8)

Hence,

$\displaystyle\mathfrak{B}(A,B;\mathbb{X})=\Gamma({A},\mathbb{X})\,\Gamma({B},\mathbb{X})\,\Gamma^{-1}({A+B},\mathbb{X}).$
(2.9)

It is obvious that $\mathbb{X}=O$ gives the original beta matrix function
[10]. Let $A$, $B$, $C$ be positive stable matrices in $\mathbb{C}^{r\times
r}$ such that $C+kI$ is invertible for all integers $k\geq 0$. Then the Gauss
hypergeometric matrix function is defined by [11]

$\displaystyle{}_{2}F_{1}(A,B;C;z)=\sum_{n=0}^{\infty}(A)_{n}(B)_{n}(C)_{n}^{-1}\frac{z^{n}}{n!}.$
(2.10)

The series (2.10) converges absolutely for $|z|<1$ and for $z=1$, if
$\alpha(A)+\alpha(B)<\beta(C)$, where $\alpha(A)=\max\\{\,\Re(z)\mid
z\in\sigma(A)\,\\}$, $\beta(A)=\min\\{\,\Re(z)\mid z\in\sigma(A)\,\\}$ and
$\beta(A)=-\alpha(-A)$.

Furthermore, if $CB=BC$ and $C$, $B$ and $C-B$ are positive stable, then for
$|z|<1$ an integral representation of (2.10) is given by [11]

$\displaystyle{{}_{2}F_{1}}(A,B;C;z)=$
$\displaystyle\left(\int_{0}^{1}(1-zt)^{-A}\,t^{B-I}(1-t)^{C-B-I}dt\right)$
$\displaystyle\times\Gamma^{-1}(B)\Gamma^{-1}(C-B)\Gamma(C).$ (2.11)

Recently, Abdalla et.al generalized the Gauss and Kummer hypergeometric matrix
function. Let $A$, $B$, $C$, $C-B$ and $\mathbb{X}$ be positive stable
matrices in $\mathbb{C}^{r\times r}$ such that $CB=BC$,
$C\mathbb{X}=\mathbb{X}C$ and $B\mathbb{X}=\mathbb{X}B$. Then the extended
Gauss hypergeometric matrix function (EGHMF) and extended Kummer
hypergeometric matrix function (EKHMF) are defined by [2]

$\displaystyle F^{(\mathbb{X})}(A,B;C;z)$ $\displaystyle=\left(\sum_{m\geq
0}(A)_{m}\,\mathfrak{B}(B+mI,C-B;\mathbb{X})\frac{z^{m}}{m!}\right)\times\Gamma(C)\Gamma^{-1}(B)\Gamma^{-1}(C-B);$
(2.12)

and

$\displaystyle\phi^{(\mathbb{X})}(B;C;z)$ $\displaystyle=\left(\sum_{m\geq
0}\,\mathfrak{B}(B+mI,C-B;\mathbb{X})\frac{z^{m}}{m!}\right)\times\Gamma(C)\Gamma^{-1}(B)\Gamma^{-1}(C-B)$
(2.13)

respectively.

Integral representation of extended matrix functions given in Equations (2.12)
and (2.13) are as follows

$\displaystyle F^{(\mathbb{X})}(A,B;C;z)$
$\displaystyle=\left(\int_{0}^{1}(1-zt)^{-A}t^{B-I}(1-t)^{C-B-I}\exp\Big{(}\frac{-\mathbb{X}}{t(1-t)}\Big{)}dt\right)$
$\displaystyle\quad\times\Gamma(C)\Gamma^{-1}(B)\Gamma^{-1}(C-B);$ (2.14)
$\displaystyle\phi^{(\mathbb{X})}(B;C;z)$
$\displaystyle=\left(\int_{0}^{1}t^{B-I}(1-t)^{C-B-I}\exp\Big{(}\frac{-\mathbb{X}}{t(1-t)}\Big{)}dt\right)\times\Gamma(C)\Gamma^{-1}(B)\Gamma^{-1}(C-B).$
(2.15)

##  3 The extended hypergeometric matrix functions of two and three variables

In this section, we introduce extended Appell matrix functions (EAMF’s) and
extended Lauricella’s matrix function (ELMF) of three variables. More
explicitly, we give the extended form of Appell matrix functions
$F_{1}(A,B,B^{\prime};C;z,w)$, $F_{2}(A,B,B^{\prime};C,C^{\prime};z,w)$ and
Lauricella matrix function of three variables
$F^{(3)}_{D}(A,B,B^{\prime},B^{\prime\prime};C;z,w,v)$, [5, 6, 7], in terms of
the extended beta matrix function. We also give here the integral
representations for these extended hypergeometric matrix functions.

Let $A$, $B$, $B^{\prime}$, $C$, $C-A$ and $\mathbb{X}$ be positive stable
matrices in $\mathbb{C}^{r\times r}$ such that $A$, $C$, $\mathbb{X}$ commutes
with each other and $CB=BC$, $CB^{\prime}=B^{\prime}C$. Then, we define
extended Appell hypergeometric matrix function
$F_{1}(A,B,B^{\prime};C;z,w;\mathbb{X})$ as

$\displaystyle F_{1}(A,B,B^{\prime};C;z,w;\mathbb{X})$
$\displaystyle=\Gamma\left(\begin{array}[]{c}C\\\
A,C-A\end{array}\right)\sum_{m,n\geq
0}\mathfrak{B}(A+(m+n)I,C-A;\mathbb{X})\,(B)_{m}(B^{\prime})_{n}\frac{z^{m}w^{n}}{m!n!},$
(3.3)

where $\Gamma\left(\begin{array}[]{c}C\\\
A,C-A\end{array}\right)=\Gamma(C)\Gamma^{-1}(A)\Gamma^{-1}(C-A)$.

For positive stable matrices $A$, $B$, $B^{\prime}$, $C$, $C^{\prime}$, $C-B$,
$C^{\prime}-B^{\prime}$ and $\mathbb{X}$ in $\mathbb{C}^{r\times r}$ such that
$B$, $B^{\prime}$, $C$, $C^{\prime}$ and $\mathbb{X}$ commutes with each
other, we define the extended Appell hypergeometric matrix function
$F_{2}(A,B,B^{\prime};C,C^{\prime};z,w;\mathbb{X})$ as

$\displaystyle F_{2}(A,B,B^{\prime};C,C^{\prime};z,w;\mathbb{X})$
$\displaystyle=\sum_{m,n\geq
0}(A)_{m+n}\mathfrak{B}(B+mI,C-B;\mathbb{X})\mathfrak{B}(B^{\prime}+nI,C^{\prime}-B^{\prime};\mathbb{X})\frac{z^{m}w^{n}}{m!n!}$
$\displaystyle\quad\times\Gamma\left(\begin{array}[]{c}C,C^{\prime}\\\
B,B^{\prime},C-B,C^{\prime}-B^{\prime}\end{array}\right).$ (3.6)

Suppose $A$, $B$, $B^{\prime}$, $B^{\prime\prime}$, $C$, $C-A$ and
$\mathbb{X}$ be positive stable matrices in $\mathbb{C}^{r\times r}$ such that
$A$, $C$, $\mathbb{X}$ commutes with each other and $CB=BC$,
$CB^{\prime}=B^{\prime}C$, $CB^{\prime\prime}=B^{\prime\prime}C$. Then, we
define the extended Lauricella hypergeometric matrix function
$F^{(3)}_{D,\mathbb{X}}(A,B,B^{\prime},B^{\prime\prime};C;z,w,v)$ as

$\displaystyle
F^{(3)}_{D,\mathbb{X}}(A,B,B^{\prime},B^{\prime\prime};C;z,w,v)$
$\displaystyle=\Gamma\left(\begin{array}[]{c}C\\\
A,C-A\end{array}\right)\sum_{m,n,p\geq
0}\mathfrak{B}(A+(m+n+p)I,C-A;\mathbb{X})(B)_{m}(B^{\prime})_{n}(B^{\prime\prime})_{p}\frac{z^{m}w^{n}v^{p}}{m!n!p!}.$
(3.9)

We now turn our attention in finding the integral representations of extended
Appell matrix functions (EAMF’s) and extended Lauricella matrix function
(ELMF) of three variables. We start with the integral representation of
$F_{1}(A,B,B^{\prime};C;z,w;\mathbb{X})$ determined in the next theorem.

######  Theorem 3.1.

Let $A$, $B$, $B^{\prime}$, $C$, $C-A$ and $\mathbb{X}$ be positive stable
matrices in $\mathbb{C}^{r\times r}$ such that $A$, $C$, $\mathbb{X}$ commutes
with each other and $CB=BC$, $CB^{\prime}=B^{\prime}C$. Then the EAMF
$F_{1}(A,B,B^{\prime};C;z,w;\mathbb{X})$ can be presented in the integral form
as

$\displaystyle F_{1}(A,B,B^{\prime};C;z,w;\mathbb{X})$
$\displaystyle=\Gamma\left(\begin{array}[]{c}C\\\
A,C-A\end{array}\right)\int_{0}^{1}u^{A-I}(1-u)^{C-A-I}(1-zu)^{-B}$ (3.12)
$\displaystyle\quad\times(1-wu)^{-B^{\prime}}\exp\Big{(}\frac{-\mathbb{X}}{u(1-u)}\Big{)}du.$
(3.13)

###### Proof.

Using the integral representation of extended beta matrix function from (2.8)
in the definition of the EAMF $F_{1}(A,B,B^{\prime};C;z,w;\mathbb{X})$, we get

$\displaystyle F_{1}(A,B,B^{\prime};C;z,w;\mathbb{X})$
$\displaystyle=\Gamma\left(\begin{array}[]{c}C\\\
A,C-A\end{array}\right)\sum_{m,n\geq
0}\int_{0}^{1}u^{A-I}(1-u)^{C-A-I}\exp\left(\frac{-\mathbb{X}}{u(1-u)}\right)$
(3.16)
$\displaystyle\quad\times(B)_{m}(B^{\prime})_{n}\frac{{(zu)}^{m}{(wu)}^{n}}{m!n!}du.$
(3.17)

Using the process discussed in [5] we can show that the sequence of matrix
functions in (3.17) is integrable and by dominated convergence theorem [8],
the summation and the integral can be interchanged in (3.17). Now applying the
matrix identity,

$\displaystyle(1-z)^{-A}=\sum_{n=0}^{\infty}(A)_{n}\frac{z^{n}}{n!},$ (3.18)

we can rewrite the Equation (3.17) as follows

$\displaystyle F_{1}(A,B,B^{\prime};C;z,w;\mathbb{X})$
$\displaystyle=\Gamma\left(\begin{array}[]{c}C\\\
A,C-A\end{array}\right)\int_{0}^{1}u^{A-I}(1-u)^{C-A-I}\exp\Big{(}\frac{-\mathbb{X}}{u(1-u)}\Big{)}$
(3.21) $\displaystyle\quad\times(1-zu)^{-B}(1-wu)^{-B^{\prime}}du.$ (3.22)

This completes the proof of Theorem 3.1. ∎

######  Theorem 3.2.

Let $A$, $B$, $B^{\prime}$, $C$, $C^{\prime}$, $C-B$, $C^{\prime}-B^{\prime}$
and $\mathbb{X}$ be positive stable matrices in $\mathbb{C}^{r\times r}$ such
that $B$, $B^{\prime}$, $C$, $C^{\prime}$ and $\mathbb{X}$ commutes with each
other. Then the EAMF $F_{2}(A,B,B^{\prime};C,C^{\prime};z,w;\mathbb{X})$
defined in (3.6) has following integral representation:

$\displaystyle F_{2}(A,B,B^{\prime};C,C^{\prime};z,w;\mathbb{X})$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}(1-zu-
wv)^{-A}u^{B-I}(1-u)^{C-B-I}v^{B^{\prime}-I}(1-v)^{C^{\prime}-B^{\prime}-I}$
$\displaystyle\quad\times\exp\left(\frac{-\mathbb{X}}{u(1-u)}-\frac{\mathbb{X}}{v(1-v)}\right)du\,dv\,\Gamma\left(\begin{array}[]{c}C,C^{\prime}\\\
B,B^{\prime},C-B,C^{\prime}-B^{\prime}\end{array}\right).$ (3.25)

###### Proof.

Using extended beta matrix function and the EAMF
$F_{2}(A,B,B^{\prime};C,C^{\prime};z,w;\mathbb{X})$ defined in (2.8) and (3.6)
respectively, we have

$\displaystyle F_{2}(A,B,B^{\prime};C,C^{\prime};z,w;\mathbb{X})$
$\displaystyle=\sum_{m,n\geq
0}\int_{0}^{1}\int_{0}^{1}(A)_{m+n}\frac{(zu)^{m}(wv)^{n}}{m!n!}u^{B-I}(1-u)^{C-B-I}v^{B^{\prime}-I}(1-v)^{C^{\prime}-B^{\prime}-I}$
$\displaystyle\quad\times\exp\left(\frac{-\mathbb{X}}{u(1-u)}-\frac{\mathbb{X}}{v(1-v)}\right)\,du\,dv\,\Gamma\left(\begin{array}[]{c}C,C^{\prime}\\\
B,B^{\prime},C-B,C^{\prime}-B^{\prime}\end{array}\right).$ (3.28)

Summation and integral in (3.28) can be interchanged by using the dominated
convergence theorem. Taking into account the summation formula [13]

$\displaystyle\sum_{N\geq 0}f(N)\frac{(z+w)^{N}}{N!}=\sum_{m,n\geq
0}f(m+n)\frac{z^{m}}{m!}\frac{w^{n}}{n!},$ (3.29)

we get

$\displaystyle F_{2}(A,B,B^{\prime};C,C^{\prime};z,w;\mathbb{X})$
$\displaystyle=\int_{0}^{1}\int_{0}^{1}\sum_{N\geq
0}(A)_{N}\frac{(zu+wv)^{N}}{N!}u^{B-I}(1-u)^{C-B-I}v^{B^{\prime}-I}(1-v)^{C^{\prime}-B^{\prime}-I}$
$\displaystyle\quad\times\exp\Big{(}\frac{-\mathbb{X}}{u(1-u)}-\frac{\mathbb{X}}{v(1-v)}\Big{)}dudv\
\Gamma\left(\begin{array}[]{c}C,C^{\prime}\\\
B,B^{\prime},C-B,C^{\prime}-B^{\prime}\end{array}\right).$ (3.32)

From (3.18) and (3.32), we get (3.25). ∎

######  Theorem 3.3.

Suppose $A$, $B$, $B^{\prime}$, $B^{\prime\prime}$, $C$, $C-A$ and
$\mathbb{X}$ be positive stable matrices in $\mathbb{C}^{r\times r}$ such that
$A$, $C$, $\mathbb{X}$ commutes with each other and $CB=BC$,
$CB^{\prime}=B^{\prime}C$, $CB^{\prime\prime}=B^{\prime\prime}C$. Then the
ELMF $F^{(3)}_{D,\mathbb{X}}(A,B,B^{\prime},B^{\prime\prime};C;z,w,v)$ defined
in (3.9) have the following integral representation:

$\displaystyle
F^{(3)}_{D,\mathbb{X}}(A,B,B^{\prime},B^{\prime\prime};C;z,w,v)$
$\displaystyle=\Gamma\left(\begin{array}[]{c}C\\\
A,C-A\end{array}\right)\int_{0}^{1}u^{A-I}(1-u)^{C-A-I}\exp\Big{(}\frac{-\mathbb{X}}{u(1-u)}\Big{)}$
(3.35)
$\displaystyle\quad\times(1-zu)^{-B}(1-wu)^{-B^{\prime}}(1-vu)^{-B^{\prime\prime}}du.$
(3.36)

###### Proof.

Equations (2.8) and (3.9) together yield

$\displaystyle
F^{(3)}_{D,\mathbb{X}}(A,B,B^{\prime},B^{\prime\prime};C;z,w,v)$
$\displaystyle=\Gamma\left(\begin{array}[]{c}C\\\
A,C-A\end{array}\right)\sum_{m,n,p\geq
0}\int_{0}^{1}u^{A-I}(1-u)^{C-A-I}\exp\Big{(}\frac{-\mathbb{X}}{u(1-u)}\Big{)}(B)_{m}(B^{\prime})_{n}$
(3.39)
$\displaystyle\quad\times(B^{\prime\prime})_{p}\frac{(uz)^{m}(uw)^{n}(uv)^{p}}{m!n!p!}.$
(3.40)

Now, using the matrix relation (3.18) and proceeding in the similar as in
Theorem 3.1, we get the required result (3.36). ∎

##  4 Fractional calculus of extended hypergeometric matrix function

The extended Riemann-Liouville fractional derivative of order $\mu$ is given
by [13]

$\displaystyle D_{z}^{\mu,p}\\{f(z)\\}$
$\displaystyle=\frac{1}{\Gamma{(-\mu)}}\int_{0}^{z}f(t)(z-t)^{-\mu-1}\exp\left(\frac{-pz^{2}}{t(z-t)}\right)dt,\quad\Re(\mu)<0,\Re(p)>0$
(4.1)

and for $m-1<\Re(\mu)<m$ (m= 1, 2,…)

$\displaystyle D_{z}^{\mu,p}\\{f(z)\\}=$
$\displaystyle\frac{d^{m}}{dz^{m}}D_{z}^{\mu-m,p}\\{f(z)\\}$ $\displaystyle=$
$\displaystyle\frac{d^{m}}{dz^{m}}\left(\frac{1}{\Gamma{(-\mu+m)}}\int_{0}^{z}f(t)(z-t)^{-\mu+m-1}\exp\left(\frac{-pz^{2}}{t(z-t)}\right)dt\right)$
(4.2)

where the path of integration is a line from $0$ to $z$ in the complex
t-plane. For the case $p=0$, we obtain the classical Riemann-Liouville
fractional derivative operator.

######  Definition 1.

Let $\mathbb{X}$ be a positive stable matrix in $\mathbb{C}^{r\times r}$ and
$\mu\in\mathbb{C}$ such that $\Re(\mu)<0$. Then, the extended Riemann-
Liouville fractional derivative of order $\mu$ is defined as follows

$\displaystyle
D_{z}^{\mu,\mathbb{X}}\\{f(z)\\}=\frac{1}{\Gamma{(-\mu)}}\int_{0}^{z}f(t)(z-t)^{-\mu-1}\exp\left(\frac{-\mathbb{X}z^{2}}{t(z-t)}\right)dt$
(4.3)

and for $m-1<\Re(\mu)<m$ (m= 1, 2,…)

$\displaystyle D_{z}^{\mu,\mathbb{X}}\\{f(z)\\}=$
$\displaystyle\frac{d^{m}}{dz^{m}}D_{z}^{\mu-m,\mathbb{X}}\\{f(t)\\}$
$\displaystyle=$
$\displaystyle\frac{d^{m}}{dz^{m}}\left[\frac{1}{\Gamma{(-\mu+m)}}\int_{0}^{z}f(t)(z-t)^{-\mu+m-1}\exp\left(\frac{-\mathbb{X}z^{2}}{t(z-t)}\right)dt\right],$
(4.4)

where the path of integration is a line from $0$ to $z$ in the complex
t-plane. For the case $\mathbb{X}=O$ we obtain the classical Riemann-Liouville
fractional derivative operator.

We start our investigation by calculating the extended fractional derivatives
of some elementary functions.

######  Theorem 4.1.

Let $A$ and $\mathbb{X}$ be positive stable matrices in $\mathbb{C}^{r\times
r}$ and $\mu\in\mathbb{C}$ such that $\Re(\mu)<0$. Then,

$\displaystyle D_{z}^{\mu,\mathbb{X}}\\{z^{A}\\}=z^{A-\mu
I}\times\frac{\mathfrak{B}(A+I,-\mu I;\mathbb{X})}{\Gamma(-\mu)}$ (4.5)

###### Proof.

From (2.8) and (4.3), we obtain

$\displaystyle D_{z}^{\mu,\mathbb{X}}\\{z^{A}\\}=$
$\displaystyle\frac{1}{\Gamma{(-\mu)}}\int_{0}^{z}t^{A}(z-t)^{-\mu-1}\exp\left(\frac{-\mathbb{X}z^{2}}{t(z-t)}\right)\,dt.$
(4.6)

Putting $\frac{t}{z}=u$ and $dt=zdu$ in the equation (4.6), we get

$\displaystyle D_{z}^{\mu,\mathbb{X}}\\{z^{A}\\}=$
$\displaystyle\frac{z^{A-\mu
I}}{\Gamma(-\mu)}\int_{0}^{1}u^{A}(1-u)^{(-\mu-1)I}\exp\left(\frac{-\mathbb{X}}{u(1-u)}\right)du$
$\displaystyle=$ $\displaystyle z^{A-\mu I}\ \frac{\mathfrak{B}(A+I,-\mu
I;\mathbb{X})}{\Gamma(-\mu)}.$

This completes the proof. ∎

######  Theorem 4.2.

Let $A$, $B$ and $\mathbb{X}$ be positive stable matrices in
$\mathbb{C}^{r\times r}$ and $\lambda,\mu\in\mathbb{C}$ such that
$\Re(\mu)>\Re(\lambda)>0$. Then,

$\displaystyle D_{z}^{\lambda-\mu,\mathbb{X}}\\{z^{A-I}(1-z)^{-B}\\}$
$\displaystyle=\Gamma(A)\Gamma^{-1}(A+(\mu-\lambda)I))z^{A+(\mu-\lambda-1)I}F^{(\mathbb{X})}(B,A;A+(\mu-\lambda)I;z).$
(4.7)

###### Proof.

Proceeding similarly as in Theorem 4.1, we get

$\displaystyle D_{z}^{\lambda-\mu,\mathbb{X}}\\{z^{A-I}(1-z)^{-B}\\}$
$\displaystyle=\frac{z^{A+(\mu-\lambda-1)I}}{\Gamma(\mu-\lambda)}\int_{0}^{1}u^{A-I}(1-uz)^{-B}(1-u)^{\mu-\lambda-1}$
$\displaystyle\quad\times\exp\left(\frac{-\mathbb{X}}{u(1-u)}\right)du.$ (4.8)

Using (2.14), we can write

$\displaystyle D_{z}^{\lambda-\mu,\mathbb{X}}\\{z^{A-I}(1-z)^{-B}\\}=$
$\displaystyle\frac{z^{A+(\mu-\lambda-1)I}}{\Gamma(\mu-\lambda)}\mathfrak{B}(A,(\mu-\lambda)I)F^{(\mathbb{X})}(B,A;A+(\mu-\lambda)I;z)$
$\displaystyle=$ $\displaystyle
z^{A+(\mu-\lambda-1)I}\Gamma(A)\Gamma^{-1}(A+(\mu-\lambda)I)\,F^{(\mathbb{X})}(B,A;A+(\mu-\lambda)I;z).$

Hence the proof is completed. ∎

The transform of several other matrix functions can be obtained easily by
using the extended Riemann-Liouville fractional derivative operators defined
in (4.3) and (4.4). Since the proofs are similar to Theorems 4.1 and 4.2, we
produce the results without proofs.

######  Theorem 4.3.

Let $A$, $B$, $B^{\prime}$and $C$ be matrices in $\mathbb{C}^{r\times r}$ such
that $AB=BA$, $AB^{\prime}=B^{\prime}A$, $AC=CA$; $A$, $C$ and $C-A$ are
positive stable and $\lambda,\mu\in\mathbb{C}$ such that
$\Re(\mu)>\Re(\lambda)>0$; $|az|<1,|bz|<1$. Then

$\displaystyle
D_{z}^{\lambda-\mu,\mathbb{X}}\\{z^{A-I}(1-az)^{-B}(1-bz)^{-B^{\prime}}\\}$
$\displaystyle=\Gamma(A)\Gamma^{-1}(A+(\mu-\lambda)I))z^{A+(\mu-\lambda-1)I}F_{1}(A,B,B^{\prime};A+(\mu-\lambda)I;az,bz;\mathbb{X}).$
(4.9)

######  Theorem 4.4.

Let $A$, $B$, $B^{\prime}$, $B^{\prime\prime}$ and $C$ be matrices in
$\mathbb{C}^{r\times r}$ such that $AC=CA$, $A$ and $C$ commute with all
remain matrices; $A$, $C$ and $C-A$ are positive stable and
$\lambda,\mu\in\mathbb{C}$ such that $\Re(\mu)>\Re(\lambda)>0$;
$|az|<1,|bz|<1$, $|cz|<1$. Then

$\displaystyle
D_{z}^{\lambda-\mu,\mathbb{X}}\\{z^{A-I}(1-az)^{-B}(1-bz)^{-B^{\prime}}(1-cz)^{-B^{\prime\prime}}\\}$
$\displaystyle=\Gamma(A)\Gamma^{-1}(A+(\mu-\lambda)I))z^{A+(\mu-\lambda-1)I}F^{(3)}_{D,\mathbb{X}}(A,B,B^{\prime},B^{\prime\prime};A+(\mu-\lambda)I;az,bz,cz).$
(4.10)

######  Theorem 4.5.

Let $A$, $B$, $B^{\prime}$ and $C$ be matrices in $\mathbb{C}^{r\times r}$
such that $B$, $B^{\prime}$, $C$ commute with each other and $B$,
$B^{\prime}$, $C$, $C-B$ and $C^{\prime}-B^{\prime}$ are positive stable. Also
let $\lambda,\mu\in\mathbb{C}$ such that $\Re(\mu)>\Re(\lambda)>0$. Then, for
$|\frac{x}{1-z}|<1$, we have

$\displaystyle
D_{z}^{\lambda-\mu,\mathbb{X}}\\{z^{B^{\prime}-I}(1-z)^{-A}F^{(\mathbb{X})}(A,B;C;\frac{x}{1-z})\\}$
$\displaystyle=\frac{z^{B^{\prime}+(\mu-\lambda-1)I}}{\Gamma(\mu-\lambda)}F_{2}(A,B,B^{\prime};C,\mu
I;x,z;\mathbb{X})\times\Gamma\left(\begin{array}[]{c}C\\\
B,C-B\end{array}\right).$ (4.13)

##  5 Generating functions

######  Theorem 5.1.

Let $A$, $B$ and $\mathbb{X}$ be positive stable matrices in
$\mathbb{C}^{r\times r}$ and $\lambda,\mu\in\mathbb{C}$ such that
$\Re(\mu)>\Re(\lambda)>0$. Then, for $|z|<|1-t|$, we have

$\displaystyle\sum_{n=0}^{\infty}\frac{(A)_{n}}{n!}F^{(\mathbb{X})}(A+nI,B;B+(\mu-\lambda)I;z)\,t^{n}=(1-t)^{-A}F^{(\mathbb{X})}{(A,B;B+(\mu-\lambda)I;\frac{z}{1-t})}.$
(5.1)

###### Proof.

Consider the matrix identity

$\displaystyle[(1-z)-t]^{-A}=(1-t)^{-A}\Big{[}1-\frac{z}{1-t}\Big{]}^{-A},$

which can be written as

$\displaystyle\sum_{n=0}^{\infty}\frac{(A)_{n}}{n!}(1-z)^{-A}\left(\frac{t}{1-z}\right)^{n}=(1-t)^{-A}\Big{[}1-\frac{z}{1-t}\Big{]}^{-A}.$
(5.2)

Now, multiplying by $z^{B-I}$ in Equation (5.2) and applying the extended
fractional derivative operator $D_{z}^{\lambda-\mu,\mathbb{X}}$, we can write

$\displaystyle\sum_{n=0}^{\infty}\frac{(A)_{n}}{n!}D_{z}^{\lambda-\mu,\mathbb{X}}\\{z^{B-I}(1-z)^{-A-nI}\\}t^{n}=(1-t)^{-A}D_{z}^{\lambda-\mu,\mathbb{X}}\\{z^{B-I}\Big{[}1-\frac{z}{1-t}\Big{]}^{-A}\\}.$

Using Theorem 4.2, we get the desired result. ∎

######  Theorem 5.2.

Let $A$, $B$, $B^{\prime}$ and $\mathbb{X}$ be positive stable matrices in
$\mathbb{C}^{r\times r}$ and $\lambda,\mu\in\mathbb{C}$ such that
$\Re(\mu)>\Re(\lambda)>0$. Then, for $|t|<\frac{1}{1+|z|}$, we have the
following generating relation

$\displaystyle\sum_{n=0}^{\infty}\frac{(B^{\prime})_{n}}{n!}F^{(\mathbb{X})}(B-nI,A;A+(\mu-\lambda)I;z)\,t^{n}$
$\displaystyle=(1-t)^{-B^{\prime}}F_{1}{(A,B,B^{\prime};A+(\mu-\lambda)I;z,\frac{-zt}{1-t};\mathbb{X})}.$

###### Proof.

Using (5.2), we have

$\displaystyle\sum_{n=0}^{\infty}\frac{(B^{\prime})_{n}}{n!}(1-z)^{n}t^{n}=(1-t)^{-B^{\prime}}\Big{[}1-\frac{-zt}{1-t}\Big{]}^{-B^{\prime}}.$
(5.3)

Now multiplying both sides of the Equation (5.3) by $z^{A-I}(1-z)^{-B}$ and
applying the extended fractional derivative operator
$D_{z}^{\lambda-\mu,\mathbb{X}}$, we get

$\displaystyle\sum_{n=0}^{\infty}\frac{(B^{\prime})_{n}}{n!}D_{z}^{\lambda-\mu,\mathbb{X}}\\{z^{A-I}(1-z)^{-(B-nI)}\\}t^{n}$
$\displaystyle=(1-t)^{-B^{\prime}}D_{z}^{\lambda-\mu,\mathbb{X}}\\{z^{A-I}(1-z)^{-B}\Big{[}1-\frac{-zt}{1-t}\Big{]}^{-B^{\prime}}\\}.$

Using the Theorems 4.2 and 4.3, we get the desired result. ∎

##  6 Conclusion

In this paper, we studied the generalized version of some Appell matrix
function and Lauricella matrix function. We obtained some generating relations
in terms of matrix functions using the extended Riemann-Liouville fractional
derivative of order $\mu$. The particular case of our results, i.e. if we take
matrices from $\mathbb{C}^{1\times 1}$, coincides with the results obtained in
[3].

## References

  * [1] M. Abdalla, A. Bakhet, Extension of Beta matrix function, Asian J. Math. Comput. Res. 9 (2016), 253–264.
  * [2] M. Abdalla, A. Bakhet, Extended Gauss hypergeometric matrix functions, Iran. J. Sci. Technol. Trans. A Sci. 42 (2018), no. 3, 1465–1470.
  * [3] M. A. Özarslan, E. Özergin, Some generating relations for extended hypergeometric functions via generalized fractional derivative operator, Math. Comput. Modelling, 52 (2010), no. 9–10, 1825–1833.
  * [4] M. A. Chaudhry, A. Qadir, H. M. Srivastava, R. B. Paris, Extended hypergeometric and confluent hypergeometric functions, Appl. Math. Comput. 159 (2004), no. 2, 589–602.
  * [5] R. Dwivedi, V. Sahai, On the hypergeometric matrix functions of two variables, Linear Multilinear Algebra 66 (2018), no. 9, 1819–1837.
  * [6] R. Dwivedi, V. Sahai, On the hypergeometric matrix functions of several variables, J. Math. Phys. 59 (2018), no. 2, 023505, 15pp.
  * [7] R. Dwivedi, V. Sahai, A note on the Appell matrix functions. Quaest. Math. 43 (2020), No. 3, pp. 321-334.
  * [8] G. B. Folland, Fourier Analysis and Its Applications, Wadsworth and Brooks, pacific Grove, CA 1992.
  * [9] J. Greene, Hypergeometric functions over finite fields, Trans. Amer. Math. Soc. 301 (1987), no. 1, 77–101.
  * [10] L. Jódar, J.C. Cortés, Some properties of gamma and beta matrix functions. Appl. Math. Lett. 11 (1998), no. 1, 89–93.
  * [11] L. Jódar, J.C. Cortés, On the hypergeometric matrix function, Proceedings of the VIIIth Symposium on Orthogonal Polynomials and Their Applications (Seville, 1997). J. Comput. Appl. Math. 99 (1998), no. 1-2, 205–217.
  * [12] A. M. Mathai, Jacobians of matrix transformations and functions of matrix argument, World Scientific Publishing Co., Inc., River Edge, NJ, 1997. xii+435 pp.
  * [13] H.M. Srivastava, H. L. Manocha, A Treatise on Generating Functions, Halsted Press (Ellis Horwood Limited, Chichester), John Wiley and Sons, New York, Chichester, Brisbane and Toronto, 1984.
  * [14] J. A. Tirao: The matrix-valued hypergeometric equation, Proc. Natl. Acad. Sci. USA, 100 (2003), 8138 –8141.
  * [15] A. Verma, On the incomplete first Appell hypergeometric matrix functions $\gamma_{1}$ and $\Gamma_{1}$, Ramanujan J., Communicated (2019).
  * [16] A. Verma, On the incomplete fourth Appell hypergeometric matrix functions $\gamma_{4}$ and $\Gamma_{4}$, Indian J. pure appl. Math , Communicated (2020).
  * [17] A. Verma, On the incomplete Srivastava’s triple hypergeometric matrix functions, Quaest. Math., Accepted (2020) DOI: 10.2989/16073606.2020.1753123.
  * [18] A. Verma, Some results on the Srivastava’s triple hypergeometric matrix functions, Asian-Eur. J. Math., Accepted (2020) DOI: 10.1142/S179355712150056X.
  * [19] A. Verma, S. Yadav, On the incomplete second Appell hypergeometric matrix functions, Linear Multilinear Algebra, Accepted (2019) DOI:10.1080/03081087.2019.1640178.