# SiO masers and 7 mm continuum emission in binary AGB stars Mira and R Aqr.

W. D. Cotton  National Radio Astronomy Observatory, 520 Edgemont Road,
Charlottesville, VA 22903, USA  E. Humphreys  M. Wittkowski  European Southern
Observatory, Karl-Schwarzschild-Str. 2, 85748 Garching bei Munchen, Germany
A. Baudry  Laboratoire d’astrophysique de Bordeaux, Université de Bordeaux,
33615 Pessac, France  A. M. S. Richards  JBCA, School of Physics and
Astronomy, University of Manchester, M13 9PL, UK  W. Vlemmings  T. Khouri
Department of Space, Earth and Environment, Chalmers University of Technology,
Onsala Space Observatory, 439 92 Onsala, Sweden  S. Etoka  JBCA, School of
Physics and Astronomy, University of Manchester, M13 9PL, UK

###### Abstract

Interactions between AGB stars and a secondary in a close orbit are one
possible explanation of why some AGB stars develop into aspherical planetary
nebulae. This study uses millimeter observations of the continuum and SiO
maser emission in several symbiotic Miras looking for evidence of an
interaction between the two stars. New JVLA observations at $\sim$45 mas
resolution are analyzed, imaging continuum and SiO masers. Two of the three
targets were detected and accurately registered continuum and line images were
derived. No clear evidence of an interaction was found between components B
and A in Mira. R Aqr has a well known jet arising from the secondary star. The
jet may be disturbing the circumstellar envelop of the AGB star or possibly
just nearly aligned with it.

###### keywords:

Stars: AGB and post-AGB – Stars:binaries:symbiotic – Radio lines:stars – Radio
continuum:stars

##  1 Introduction

A substantial fraction of main sequence stars, up to 50%, are known to exist
in systems with multiple components. The importance of close binaries on
stellar evolution has been recognized lately, but their effects on the close
circumstellar environment is under active study. The influence of a companion
on the region near to the surface of Asymptotic Giant Branch (AGB) stars may
play a pivotal role in the formation of aspherical Planetary Nebulae (PNe).
Observations of molecules in the winds of a sample of AGB stars (Ramstedt et
al., 2014, 2017, 2018; Doan et al., 2017; Decin et al., 2020) show a variety
of non spherical shapes which are attributed to the interactions with
companion stars. Similar shaping might affect the episodes of mass loss
leading to the formation of PNe.

SiO masers are usually found within the pulsating atmospheres of AGB stars at
about 2 stellar radii, at about the distance at which silicate and other dust
forms in the outflows (Karovicova et al., 2013; Gobrecht et al., 2016; Höfner
& Freytag, 2019). These masers generally appear in the form of a ring centered
on the photosphere and reveal motions in the inner envelope (Boboltz et al.,
1997; Hollis et al., 1997b; Cotton et al., 2004, 2006; Ragland et al., 2008).
Single dish monitoring of semi-regular and Mira AGB stars by Gómez-Garrido et
al. (2020) show large and rapid variations of the SiO masers in the semi-
regular giants RX Boo and RT Vir and slower variations in several Miras.

One of the current unknowns in stellar evolution is the shaping mechanism from
asymptotic giant branch (AGB) stars to planetary nebulae. On the one hand,
both recent optical spectro–polarimetry and radio/submm maser observations
indicate unexpectedly high magnetic fields near or at the stellar surface
(Lèbre et al., 2014; Vlemmings, 2018) and references therein. On the other
hand, ALMA observations indicate the presence of companions around stars that
were previously believed to be isolated e.g. the iconic ALMA observations of R
Sculptoris; (Maercker et al., 2012). Whether binarity or magnetic fields
dominate the evolutionary shaping mechanism is a matter of current debate.

One of the best places to study how the out-flowing material is shaped is at,
or near, the stellar surface. Long baseline ALMA observations of the surfaces
of AGB stars Mira and W Hya reveal hot-spots that likely influence and cause
asymmetry in the near–circumstellar environment (Vlemmings et al., 2015,
2017). Towards binary system Mira AB, we have also recently discovered signs
of the influence of the white dwarf Mira B at only a few stellar radii from
AGB star Mira A (Humphreys et al., in prep). In this case, using ALMA long
baseline data we found that SiO masers trace a portion of a bubble wall formed
by the interacting Mira AB winds, also seen as a spiral plume in SiO thermal
gas. In this paper, we present observations of continuum and SiO maser lines
using the JVLA of Mira and two other AGB stars with known companions, W Aql
and R Aqr. Our aims are to study the effect of the companions on the observed
maser spot distribution.

Cotton et al. (2004, 2006) give a sequence of VLBA images of the SiO masers in
the atmosphere of Mira. These appear as partial rings with diameters in the
range 60 to 75 mas. The diameter of the $\nu$=1,J=1-0 maser ring tends to be
slightly larger than that of the $\nu$=2,J=1-0 masers. Mennesson et al. (2002)
give a uniform disk diameter of the photosphere in K band of 28.79 $\pm$ 0.10
mas. Wittkowski et al. (2016) give a diameter of 28.5 $\pm$ 1.5 mas at 2.25
$\mu$m. Perrin et al. (2020) give the photospheric diameter in H band of 21.1
mas. The distance to Mira is given by Haniff et al. (1995) as 110 $\pm$ 9 pc.
at which distance 1” corresponds to 110 AU. The Hipparcos distance (van
Leeuwen, 2007) is 92$\pm$11 pc.

W Aquilae (W Aql) is a binary S–type AGB star sometimes showing SiO maser
emission. SiO maser emission observed towards W Aql (Nakashima & Deguchi,
2007) disappeared for a time around 2011 (Ramstedt et al., 2012), which could
be due to disruption by an additional companion (the known companion lies at
about 0.46 arcseconds from the S–type AGB star, about 180 AU at a distance to
W Aql of 395 pc). Peak SiO maser emission can reach 21 Jy (single–dish;
Ramstedt et al. (2012)).

R Aquarii (R Aqr) is a symbiotic binary and a well–known host of millimeter
and submillimeter SiO masers (Boboltz et al., 1997; Hollis et al., 1997b; Gray
et al., 1998; Cotton et al., 2004, 2006; Ragland et al., 2008). Cotton et al.
(2004, 2006) give a sequence of VLBA images of the SiO masers in the
atmosphere of R Aqr with ring diameters of 31 to 33 mas. Infrared measurements
of the photospheric size show a range of values between 11.2 and 18.4 mas
(Tuthill et al., 2000; Mennesson et al., 2002; Ragland et al., 2008;
Wittkowski et al., 2016).

The maser emission is variable but peak values in single–dish observations lie
in the range up to 300 Jy. R Aqr consists of a Mira variable and an accreting,
hot companion with a remarkable jet outflow. HST observations of the jet shows
considerable change over a several year timescale (Hollis et al., 1997a). In
2014 the stars were separated by 45 mas (Schmid et al., 2017) or 9.8 AU at the
distance to R Aqr of 218 pc (Min et al., 2014). ALMA observations of
Bujarrabal et al. (2018) detect both stars with a bridge of material joining
them, likely material flowing from the AGB primary to the accretion disk
around the WD secondary. Cotton et al. (2004) show an unusual rotating ring of
masers in January 2001\. This ring was interpreted as arising in a rotating
equatorial disk of material.

##  2 Observations

The observations were made on the JVLA near 7 mm wavelength during 2019 August
12 from 06:00:00 UT to 13:00:00 UT under project code 19A-220 in the most
extended (“A”) configuration. The target stars (astrometric/bandpass
calibrators) were W Aquilae (J1939-1525), R Aquarii (J2348-1631) and Mira, AKA
omicron Ceti, (J0217+0144). The pulsation phase for W Aquilae was 0.27, for R
Aquarii was 0.21 and 0.79 for Mira. The photometric calibrator was 3C48.

The observing sequence consisted of cycling between one of the targets (7 min)
and its calibrator (2 min) for an hour and then proceeding to another target.
Reference pointing was used with a pointing measurement made on the calibrator
at the beginning of each sequence. The integration time was 2 seconds and only
the parallel hand data (RR,LL) were recorded.

The masers observed were the SiO $\nu$=2,J=1-0 and $\nu$=1,J=1-0 transitions
at 42.820587 and 43.122027 GHz respectively. Each transition was covered by
256 channels of 62.5 kHz bandwidth($\sim$0.9 km/s after Hanning). The center
LSR velocity of the observation was 47 km/sec for Mira and -22 km/sec for R
Aqr and W Aql for which the true systemic velocity is approximately -26 and
-18 km/sec respectively. The continuum observations consisted of 58 spectral
windows in the range of 41.0 to 48.7 GHz using 7424 $\times$ 1 MHz channels;
spectral regions containing the masers were not included.

##  3 Calibration

The conditions of these observations, Summertime, high frequency and high
resolution render the proper calibration of the data more difficult. In order
to allow accurate comparison of the continuum and line images they must be
aligned both photometrically and especially astrometrically. The relatively
poor phase coherence of the data make this more difficult. Photometric
calibration is complicated by the relatively poor sensitivity and strong
resolution of the photometric calibrator by the $\sim$45 mas resolution of the
given data. Fortunately, the masers in two of the targets are quite strong and
allow relative phase calibration following the general approach of Reid &
Menten (1990, 1997, 2007), adapted to the wide-band nature of the present
data.

Data were translated from the archive format to AIPS format and were processed
in the Obit package (Cotton, 2008). Hanning smoothing was used on the line
data to suppress the artifacts resulting from the interaction of the strong,
narrow maser signals and the finite delay range of the correlator.

###  3.1 Amplitude

The traditional amplitude calibration scheme for connected element
interferometers is to compare the gain solutions of the astrometric and
photometric calibrators to infer the flux density of the astrometric
calibrator. This works well in the high signal-to-noise regime with well
modeled and only marginally resolved calibrators but can produce poor results
with low SNR data with strongly resolved photometric calibrators as is the
case here. All of the centimeter wavelength stable photometric calibration
sources, including 3C48 used in this work, are strongly resolved at the
$\sim$45 mas resolution of this data.

We have adopted an alternate approach of using an initial, assumed flux
density for the astrometric calibrator and apply this calibration to the
photometric calibrator. The ratio of the apparent integrated flux density of
the derived image of the photometric calibrator to its true flux density gives
the factor needed to correct the assumed flux density of the astrometric
calibrator. Astrometric calibrators are generally physically very small
sources, hence variable, but are well modeled by a point at the current
resolution. This approach has the advantage that a super accurate structural
model of the photometric calibrator is not needed but only an accurate total
flux density. The adopted flux density of 3C48 at 40.98 GHz was 0.611 Jy and
at 48.6 GHz was 0.506 Jy (Perley & Butler, 2017). The systematic error in the
amplitude scale is estimated to be 20%.

###  3.2 Phase

The phase calibration approach of Reid & Menten (1990, 1997) is to use the
phase calibration derived from a self calibration of a strong, simple maser to
calibrate the phase of the continuum data. This can work well in the present
case of AGB stars which may have very strong, non-thermal masers but very weak
thermal emission from the photosphere.

The original scheme of Reid & Menten (1990, 1997) was appropriate for a narrow
band system for which a standard group delay calibration was adequate and did
not need to be determined from the data-set in question. The continuum
bandwidth of the current data is sufficiently wide that corrections to the
group delay need to be determined and applied to the data.

The “phase” corrections derived from the self calibration of a maser are
largely due to tropospheric refraction, hence are narrow band samples of a
group delay function and this needs to be taken into account in transferring
the maser calibration to the continuum data which extends for several GHz from
the frequency of the maser. Furthermore, the signal from a given maser spot is
very narrow band and inadequate to measure the relevant group delay. The
calibration method used is similar to that described by Matthews et al. (2015,
2018).

We have adopted the following scheme.

  1. 1.

System temperature calibration. Amplitude corrections were determined using
program SYGain from the on–line calibration measurements.

  2. 2.

Group delay calibration. The continuum observations of the phase reference
sources were used to determine the group delay error (program Calib). These
calibrations were transferred to the maser line data and used to calibrate
both data sets.

  3. 3.

Bandpass calibration. The observations of the phase reference sources were
used to derive a bandpass response function for each of the line and continuum
data for each target star using program BPass.

  4. 4.

External calibration. The amplitude calibration scheme described in Section
3.1 was used to derive the flux densities of the astrometric calibrators.
These were then used to obtain complex gain corrections (program Calib) which
were applied to the data.

  5. 5.

Self calibrate a maser. A strong, simple maser appearing in a few line
channels was imaged by SCMap for each target using self calibration. After
several cycles of phase only calibration, an amplitude and phase self
calibration was done.

  6. 6.

Transfer phase to continuum data. The time sequences of gain corrections from
the maser self calibration were transferred to the continuum data
extrapolating phase in frequency using frequency scaling. This is done with
program SnCpy.

  7. 7.

Doppler corrections. Program CVel was used to correct the spectroscopic data
for the Earth’s motion.

  8. 8.

Image. The continuum data were imaged with no further calibration producing a
wide-band image using program MFImage. Channels in which maser emission might
appear were excluded from the continuum imaging. A spectral cube was imaged
(via Imager) from the maser data for each transition using the calibration
from the self calibration.

Use of self calibration to set the coordinates of the derived images results
in a large uncertainty of the absolute positions. However, the resulting
images should have accurate relative positions of both masers and continuum
emission with respect to the reference maser.

##  4 Results

No maser or continuum emission was detected from W Aql but strong masers were
detected in Mira and R Aqr and were used to calibrate the continuum. Frequency
channels with maser emission are identified from long term, scalar averaged,
visibility spectra from individual baselines. The mean channel (0.9 km sec-1 )
visibility amplitude on a given baseline for W Aql is 18 mJy and 5 times this
should be easily detectable, a plausible upper limit to a maser peak is 90
mJy. Pimpanuwat et al. (in prep.) report that ALMA observations of the ATOMIUM
111https://fys.kuleuven.be/ster/research-projects/aerosol/atomium/atomium
project detected thermal but not maser emissions of several 1.2 mm SiO J=5-4
and J=6-5 transitions in W Aql during June and July 2019.

Many maser features are well separated in velocity and are quite small;
centroids can be measured to a small fraction of the resolution. The maser
rings can therefore be defined with an accuracy much better than the size of
the stellar photosphere. The noise contribution to the position uncertainty of
an unresolved, isolated component is $\sim\theta_{FWHM}/(2\times{\rm SNR})$
where $\theta_{FWHM}$ is the full width at half maximum of the resolution and
${\rm SNR}$ is the signal-to-noise ratio. Since the spatial and frequency
resolution of this data is such that masers overlap in space and velocity, a
list of “masers” was derived by fitting one or two Gaussians to the image in
alternate channels in each of the transitions.

###  4.1 Mira

The off–source channel noise in the maser cubes is 3.5 and 3.0 mJy beam-1 for
$\nu$=1,J=1-0 and $\nu$=2,J=1-0 respectively and the CLEAN restoring beams
were 51 $\times$ 45 mas at position angle 41∘ and 50 $\times$ 40 mas at
position angle -27∘. The integrated spectra of the SiO masers in Mira are
shown in Figure 1.

The locations of the masers are plotted on the gray scale of the continuum in
Figure 2. The velocities and relative strengths of the masers are shown in
Figure 3. The SNR of the masers seen in Figure 1 exceeds 100 so the
uncertainty due to noise of the maser locations is less than 0.25 mas, much
less than the maser ring size of 60 - 70 mas (Cotton et al., 2004) or the
photospheric size of 28.79 mas (Mennesson et al., 2002). The velocity ranges
of the masers in the two transitions isn’t the same but this is also the case
at several epochs shown in Cotton et al. (2006) and is common in the single
dish spectra of Gómez-Garrido et al. (2020). The continuum of the secondary,
Mira B, is clearly detected but no SiO was detected except at the location of
Mira A. No continuum was detected away from the two stars.

The off-source RMS in the continuum image is 28 $\mu$Jy beam-1, and the CLEAN
restoring beam is 41 $\times$ 39 mas at position angle -79∘. The peak of Mira
B is 434 $\pm$ 3 mas (53 $\pm$ 0.3 AU) from Mira A at an orientation of 95∘
from north through east. The peak flux density of Mira A is 2.7 $\pm$ 0.6 mJy.
The total flux density of Mira A is 5.0 $\pm$ 1.0 mJy with a deconvolved
Gaussian size (FWHM) of 39.4 $\pm$ 1.5 mas $\times$ 32.9 $\pm$ 1.5 mas at
position angle -56.1∘ $\pm$ 4∘.

Mira A being only moderately resolved and its profile, hence area, being not
well constrained by the present data, two possible extremes are a uniform disk
and a Gaussian. We directly use the factor of 1.6 in (Pearson, 1999) to
convert the Gaussian size to an equivalent uniform disk size. We thus get for
this extreme possibility a uniform disk diameter of 63.0$\times$52.6 mas. The
brightness temperature of Mira A calculating its area assuming this uniform
disk and a total flux density of 5.0 mJy is 1300 $\pm$ 270 K. The brightness
temperature calculating its area assuming a Gaussian profile with the flux
density of 5.0 mJy is 2320 $\pm$ 480K. The true brightness temperature is
likely between 1300 and 2300 K. Our total flux density was obtained by an
intergral over the image pixels and its error is dominated by the 20%
estimated uncertainty in the flux density scale. The errors on the Gaussian
size include the effects of correlated noise in the image (Condon, 1997). The
true uncertainty in the brightness temperature is dominated by the uncertainty
in the effective area of the source.

The peak flux density of Mira B is 0.62 $\pm$ 0.14 mJy with an integrated
value of 0.64 $\pm$ 0.14 mJy. This component is at most marginally resolved
with an upper limit of 20 mas.

Figure 1: Spectrum of Mira SiO masers, “+” is $\nu$=1,J=1-0, and “x” is
$\nu$=2,J=1-0. Figure 2: Continuum reverse gray-scale of Mira with “+” marking
the locations of $\nu$=1,J=1-0 SiO masers and “x” for $\nu$=2,J=1-0 masers. A
scalebar for the grayscale in $\mu$Jy Beam-1 is given at the top. The
resolution is shown in the box in the lower left corner. Figure 3: Color coded
velocities of SiO masers surrounding Mira A with “+” marking the locations of
$\nu$=1,J=1-0 SiO masers and “x” for $\nu$=2,J=1-0 masers. Velocities are
indicated by the color with a scale bar to the right. The size of each symbol
is proportional to the $\sqrt{\quad}$ of the peak flux density.

###  4.2 R Aquarii

The off–source channel noise in the maser cubes is 4.4 and 3.9 mJy beam-1 for
$\nu$=1,J=1-0 and $\nu$=2,J=1-0 respectively and the CLEAN restoring beams
were 59 $\times$ 41 mas at position angle -4∘. The integrated spectra of the
SiO masers in R Aqr are shown in Figure 4.

The locations of the masers are plotted on the gray scale image of the
continuum in Figure 5. The velocities and relative strengths of the masers are
shown in Figure 6. As seen from Figure 4, the SNR of the masers exceed 100 so
the uncertainty of the maser locations due to image noise is less than 0.3
mas.

The off–source RMS in the continuum image is 36 $\mu$Jy beam-1 and the CLEAN
restoring beam is 50 $\times$ 36 mas at position angle -1∘. The continuum
component in which the masers are imbedded has a peak flux density of 6.1
$\pm$ 1.2 mJy and an integrated flux density of 10.8 $\pm$ 2.2 mJy and a
deconvolved Gaussian size of 58.9 $\pm$ 1.6 mas $\times$ 16.9 $\pm$ 1.8 mas at
a position angle of 17∘ $\pm$ 1∘. The Gaussian brightness temperature of this
component is 6,550 $\pm$ 1490 K. For the equivalent uniform disk, the
brightness temperature is 3,690 $\pm$ 840 K. The integrated flux density of
the emission visible in Figure 5 is 40 mJy. No continuum or masers were
detected outside of the region shown in Figure 5.

Figure 4: Spectrum of R Aquarii SiO masers, “+” is $\nu$=1,J=1-0, and “x” is
$\nu$=2,J=1-0. Figure 5: Continuum reverse grayscale of R Aqr with “+” marking
the locations of $\nu$=1,J=1-0 SiO masers and “x” for $\nu$=2,J=1-0 masers. A
scalebar for the grayscale in $\mu$Jy Beam-1 is given at the top. Labels mark
features approximately corresponding to those in Schmid et al. (2017). The
circle shows the location of the secondary as estimated from the orbit of
Bujarrabal et al. (2018). The resolution is shown in the box in the lower left
corner Figure 6: Color coded velocities of SiO masers surrounding R Aqr with
“+” marking the locations of $\nu$=1,J=1-0 SiO masers and “x” for
$\nu$=2,J=1-0 masers. Velocities are indicated by the color with a scale bar
to the right. The size of each symbol is proportional to the $\sqrt{\quad}$ of
the peak flux density.

##  5 Discussion

###  5.1 Mira

The continuum emission from Mira A is slightly elongated 40∘ from the
direction of Mira B. Otherwise there is no hint in the present results of an
interaction. While not well defined, the maser ring is basically round.
Numerous measurements of the 7 mm SiO masers are available in the literature.
The maser rings seen at much higher resolution in Cotton et al. (2004) and
Cotton et al. (2006) are even less well defined but also show little hint of
an interaction. While there are variations in maser velocity around the ring
seen in Figure 3, this is commonly seen in AGB maser rings and shows no
obvious effect of the companion. SiO masers tracing the interaction of the AGB
star with a close companion was reported in $\pi^{1}$ Gru (Homan et al., 2020)
although with a closer companion (6 AU) than Mira B is to A.

Reid & Menten (2007) give a uniform disk size for Mira A of 54 $\pm$ 4 mas
$\times$ 50 $\pm$ 4 mas at position angle -30∘ $\pm$ 50 with a flux density of
4.8 $\pm$ 0.2 mJy and a brightness temperature of 1680 $\pm$ 250 K in 2000
October. This is somewhat smaller than our uniform disk equivalent of
63.0$\times$52.6 mas. The Reid & Menten (2007) flux density of 4.8 mJy is in
good agreement with our value of 5.0 mJy. Our larger size results in a lower
estimated uniform disk brightness temperature, 1300 $\pm$ 270 K at 45 GHz;
however, this difference is within the errors. The alignment of the SiO maser
ring with Mira A in Reid & Menten (2007) generally agrees with that in Figure
2; centered on the photosphere. They do not detect Mira B.

Observations of Mira A in 2014 February by Matthews et al. (2015) give an
elliptical Gaussian fit at 46 GHz of 37.5 $\pm$ 2.1 $\times$ 31.7 $\pm$ 2.0
mas at position angle 147∘ $\pm$ 6∘ in excellent agreement with our 2019
August results. However, the flux density quoted by Matthews et al. (2015),
7.6 $\pm$ 1.5 mJy, is 1.5 $\sigma$ higher than our 5.0 mJy and their uniform
disk brightness temperature of 2110 $\pm$ 440 K is correspondingly higher than
our 1300 $\pm$ 270 K. The difference between the current results and those of
Matthews et al. (2015) is 1.5 $\sigma$. Matthews et al. (2018) presents
evidence for variations in the sizes and shapes of several other Mira
variables so the differences between our results and those of Reid & Menten
(2007) and Matthews et al. (2015) may simply reflect changes in Mira A.

The observations of Matthews et al. (2015) detected Mira B and give a
deconvolved Gaussian size at 46 GHz of 18 $\pm$ 5 $\times$ 1 $\pm$ 6 mas at
position angle 172∘ $\pm$ 30∘ with a flux density of 0.97 $\pm$ 0.2 mJy. This
is in good agreement with our upper limit of 20 mas although the star is at
most marginally resolved in both cases. Since Mira B is thought to be a white
dwarf, any resolution of this object is very unlikely to be of the star itself
but more likely an accretion disk around it.

While we detect no emission away from Mira A and B, Wong et al. (2016) show
ALMA observations of a plume of molecular gas in an arc from Mira A up to 3”
in the direction away from Mira B but do not associate it with an effect by
Mira B. They also report that the SiO density drops sharply past 4 stellar
radii.

ALMA observations by Kamiński et al. (2016) of gas phase AlO, the precursor
molecule to Al2O3 in aluminum dust, shows that the density peaks at about
2R$\star$. This is near the size of the radio photosphere and the SiO maser
ring seen in Figure 2.

###  5.2 R Aquarii

The orbit of the secondary star in R Aqr has been the subject of debate for a
number of years. Hollis et al. (1997b) suggest that an extended feature seen
in a VLA 7 mm image corresponded to the secondary. However, Schmid et al.
(2017) dispute this interpretation and on the basis of high spatial resolution
Hα imaging suggest another. A spectroscopic orbit has been determined by
Gromadzki & Mikołajewska (2009) who derive an eccentric orbit ($e$=0.25) with
a 43.6 year period but cannot determine where on the sky the secondary
appears. Detection of the secondary star by Bujarrabal et al. (2018) allowed
them to update the orbit of Gromadzki & Mikołajewska (2009) giving an estimate
of the location at the epoch of the observations presented here. The region of
emission enclosing the stars is shown in Figure 7 with the estimated location
of the white dwarf shown as a circle. As for Mira, Figure 6 shows no clear
signature of an effect of the secondary on the velocity structure of the ring.

However, ALMA observations by Bujarrabal et al. (2018) show plumes of CO gas
which they intreprete as arising from the outflow from the AGB star being
focused into its orbital plane by the white dwarf companion maximizing the
accretion onto the secondary. They also detect continuum emission between the
stars with a spectral index suggesting dust. Such structure would be obscured
at 7 mm wavelength in 2019 by the inner jet emission seen in Figure 7.

The continuum emission seen in the ALMA 0.8 mm observations are not expected
to be the same as the 7 mm continuum emission as the former is dominated by
the photosphere and the dust forming around the star and the latter is from
free-free emission from ionized material in the circumstellar envelop.

Figure 7: Like Figure 5 but a closeup of the emission enclosing the stars and
a scalebar labeled in mJy beam-1 The circle shows the location of the
secondary as estimated from the orbit of Bujarrabal et al. (2018).

This system contains a very extended and complex jet to the north–east and
south–west which is presumed to arise from the secondary star (Hollis et al.,
1997a, b; Schmid et al., 2017). The jet material is quite visible in the Hα
images presented in Schmid et al. (2017) which also include the potential
secondary. The 7 mm image in Figure 5 is strikingly similar to several of the
Hα features south-west of the AGB star shown in Schmid et al. (2017),
especially their features C, ASW, BSW and CSW and are approximately indicated
in Figure 5. These features in Figure 5, with the partial exception of C, are
presumably related to the jet. This correspondance is puzzling given the 5
years between the observations and the very short recombination times of these
features estimated by Schmid et al. (2017). This suggests either longer
recombination times or a continuing source of ionization.

The continuum feature associated with the SiO masers in Figures 5 and 7 is
elongated in the general direction of the jet and is not centered on the maser
ring but rather on the companion white dwarf. This suggests that this feature
also encloses the secondary and some of its emission is from the jet. The
Gaussian brightness temperature of this feature, 6,550 $\pm$ 970 K is greatly
in excess of that expected for the photosphere of an AGB star further
suggesting that much of the emission is from the accretion disk/jet which is
nearly aligned on the sky with the AGB star.

The circumstellar masers seen in Figures 5 and 7 show little evidence of being
affected by the jet or the secondary star. The VLBI image of the masers in R
Aqr shown in Ragland et al. (2008) show masers streaming off in the direction
of the north–east jet in 2006 September but those in Cotton et al. (2004,
2006) do not.

The continuum emission see in Figures 5 and 7 is similar in shape and
orientation to that reported by Hollis et al. (1997b) from 1996, a half
orbital period earlier, but with a different offset from the SiO masers
marking the atmosphere of the AGB star. This suggests that in both cases, the
emission is dominated by the jet from the accretion disk around the white
dwarf and that the image of Hollis et al. (1997b) does show the location of
the secondary.

Unfortunately, the poor absolute astrometry resulting from the atmospheric
instability in the current data prevent them from being very useful
quantitatively in refining the orbit.

The distributions of masers in AGB atmospheres are extremely variable as they
track conditions suitable for maser amplification rather than following
concentrations of material. Observations several months apart, as shown in the
references for VLBI observations usually show very different distributions.
The frequent monitoring of the Mira TX Cam (Gonidakis et al., 2013) shows very
turbulent behavior.

The nature of the SiO maser distribution in January 2001 from Cotton et al.
(2004) was of a very different character than at much later epochs in showing
a high degree of symmetry across the ring. In April 2001 the general velocity
pattern of the masers was consistent with a general rotation but without much
evidence for an equatorial disk. Hollis et al. (2001) claim some evidence for
differential rotation in Dec 2000 as well as earlier with a pole of the
rotation at a position angle of 150∘. By 2006 (Ragland et al., 2008) the
velocity structure was incompatible with the differential rotation of
2000-2001 and the 2019 velocity structure in Figure 6 is also incompatible.
The 2000-2001 event may have been a transient feature although the orientation
of the apparent pole is difficult to explain by an interaction with either the
gravity of, or the jet coming from, the secondary. The jet is closer to the
direction of the pole of the putative rotation rather than orthogonal to it
and the orbit of Bujarrabal et al. (2018) shows a substantial separation
between the two stars in 2001.

##  6 Conclusions

Observations of 7 mm continuum and SiO masers in three symbotic AGB stars are
reported. The continuum data were phase referenced to the masers allowing
sensitive imaging and accurate registration. Techniques for photometric
calibration and this phase referencing were developed to exploit these data.

No masers were detected in W Aquilae which prevented a sensitive continuum
image from being made. Both Mira and R Aquarii have strong masers which allows
deep continuum imaging. The image of Mira shows the masers well centered on
the photosphere of Mira A and the white dwarf Mira B was well detected. No
clear evidence of an interaction between the white dwarf and the envelope of
Mira A was found.

The inner portion of the jet in R Aquarii is visible in the continuum image
and the radio emission around the AGB star is extended in the direction of the
jet. The continuum emission is centered on the location of the white dwarf
estimated from the orbit of Bujarrabal et al. (2018) rather than on the AGB
star as indicated by the SiO masers. That, plus the high brightness
temperature of this feature, indicates that it is dominated by emission from
the accretion disk + jets associated with the white dwarf. Bujarrabal et al.
(2018) also shows gas apparently flowing from the AGB star onto the secondary.
The alignment of the continuum peak of R Aqr with the predicted location from
the orbit of the secondary by Bujarrabal et al. (2018) appears to support
their orbit.

###### Acknowledgements.

We thank B. Pimpanuwat/ATOMIUM consortium for the J=5-4 and J=6-5
measurements. This work made use of the AAVSO International Database
contributed by observers worldwide.

## References

  * Boboltz et al. (1997) Boboltz, D. A., Diamond, P. J., & Kemball, A. J. 1997, ApJ, 487, L147
  * Bujarrabal et al. (2018) Bujarrabal, V., Alcolea, J., Mikołajewska, J., Castro-Carrizo, A., & Ramstedt, S. 2018, A&A, 616, L3
  * Condon (1997) Condon, J. J. 1997, PASP, 109, 166
  * Cotton (2008) Cotton, W. D. 2008, PASP, 120, 439
  * Cotton et al. (2004) Cotton, W. D., Mennesson, B., Diamond, P. J., et al. 2004, A&A, 414, 275
  * Cotton et al. (2006) Cotton, W. D., Vlemmings, W., Mennesson, B., et al. 2006, A&A, 456, 339
  * Decin et al. (2020) Decin, L., Montargès, M., Richards, A. M. S., et al. 2020, Science, 369, 1497
  * Doan et al. (2017) Doan, L., Ramstedt, S., Vlemmings, W. H. T., et al. 2017, A&A, 605, A28
  * Gobrecht et al. (2016) Gobrecht, D., Cherchneff, I., Sarangi, A., Plane, J. M. C., & Bromley, S. T. 2016, A&A, 585, A6
  * Gómez-Garrido et al. (2020) Gómez-Garrido, M., Bujarrabal, V., Alcolea, J., et al. 2020, A&A, 642, A213
  * Gonidakis et al. (2013) Gonidakis, I., Diamond, P. J., & Kemball, A. J. 2013, MNRAS, 433, 3133
  * Gray et al. (1998) Gray, M. D., Ivison, R. J., Humphreys, E. M. L., & Yates, J. A. 1998, MNRAS, 295, 970
  * Gromadzki & Mikołajewska (2009) Gromadzki, M., & Mikołajewska, J. 2009, A&A, 495, 931
  * Haniff et al. (1995) Haniff, C. A., Scholz, M., & Tuthill, P. G. 1995, MNRAS, 276, 640
  * Höfner & Freytag (2019) Höfner, S., & Freytag, B. 2019, A&A, 623, A158
  * Hollis et al. (2001) Hollis, J. M., Boboltz, D. A., Pedelty, J. A., White, S. M., & Forster, J. R. 2001, ApJ, 559, L37
  * Hollis et al. (1997a) Hollis, J. M., Lyon, R. G., Dorband, J. E., & Feibelman, W. A. 1997a, ApJ, 475, 231
  * Hollis et al. (1997b) Hollis, J. M., Pedelty, J. A., & Lyon, R. G. 1997b, ApJ, 482, L85
  * Homan et al. (2020) Homan, W., Montargès, M., Pimpanuwat, B., et al. 2020, A&A, 644, A61
  * Kamiński et al. (2016) Kamiński, T., Wong, K. T., Schmidt, M. R., et al. 2016, A&A, 592, A42
  * Karovicova et al. (2013) Karovicova, I., Wittkowski, M., Ohnaka, K., et al. 2013, A&A, 560, A75
  * Lèbre et al. (2014) Lèbre, A., Aurière, M., Fabas, N., et al. 2014, A&A, 561, A85
  * Maercker et al. (2012) Maercker, M., Mohamed, S., Vlemmings, W. H. T., et al. 2012, Nature, 490, 232
  * Matthews et al. (2015) Matthews, L. D., Reid, M. J., & Menten, K. M. 2015, ApJ, 808, 36
  * Matthews et al. (2018) Matthews, L. D., Reid, M. J., Menten, K. M., & Akiyama, K. 2018, AJ, 156, 15
  * Mennesson et al. (2002) Mennesson, B., Perrin, G., Chagnon, G., et al. 2002, ApJ, 579, 446
  * Min et al. (2014) Min, C., Matsumoto, N., Kim, M. K., et al. 2014, PASJ, 66, 38
  * Nakashima & Deguchi (2007) Nakashima, J.-i., & Deguchi, S. 2007, ApJ, 669, 446
  * Pearson (1999) Pearson, T. J. 1999, in Astronomical Society of the Pacific Conference Series, Vol. 180, Synthesis Imaging in Radio Astronomy II, ed. G. B. Taylor, C. L. Carilli, & R. A. Perley, 335
  * Perley & Butler (2017) Perley, R. A., & Butler, B. J. 2017, ApJS, 230, 7
  * Perrin et al. (2020) Perrin, G., Ridgway, S. T., Lacour, S., et al. 2020, A&A, 642, A82
  * Ragland et al. (2008) Ragland, S., Le Coroller, H., Pluzhnik, E., et al. 2008, ApJ, 679, 746
  * Ramstedt et al. (2018) Ramstedt, S., Mohamed, S., Olander, T., et al. 2018, A&A, 616, A61
  * Ramstedt et al. (2012) Ramstedt, S., Vlemmings, W., Mohamed, S., Choi, Y. K., & Olofsson, H. 2012, in Cosmic Masers - from OH to H0, ed. R. S. Booth, W. H. T. Vlemmings, & E. M. L. Humphreys, Vol. 287, 260–261
  * Ramstedt et al. (2014) Ramstedt, S., Mohamed, S., Vlemmings, W. H. T., et al. 2014, A&A, 570, L14
  * Ramstedt et al. (2017) —. 2017, A&A, 605, A126
  * Reid & Menten (1990) Reid, M. J., & Menten, K. M. 1990, ApJ, 370, L51
  * Reid & Menten (1997) —. 1997, ApJ, 476, 327
  * Reid & Menten (2007) —. 2007, ApJ, 671, 2068
  * Schmid et al. (2017) Schmid, H. M., Bazzon, A., Milli, J., et al. 2017, A&A, 602, A53
  * Tuthill et al. (2000) Tuthill, P. G., Danchi, W. C., Hale, D. S., Monnier, J. D., & Townes, C. H. 2000, ApJ, 534, 907
  * van Leeuwen (2007) van Leeuwen, F. 2007, A&A, 474, 653
  * Vlemmings et al. (2017) Vlemmings, W., Khouri, T., O’Gorman, E., et al. 2017, Nature Astronomy, 1, 848
  * Vlemmings (2018) Vlemmings, W. H. T. 2018, Contributions of the Astronomical Observatory Skalnate Pleso, 48, 187
  * Vlemmings et al. (2015) Vlemmings, W. H. T., Ramstedt, S., O’Gorman, E., et al. 2015, A&A, 577, L4
  * Wittkowski et al. (2016) Wittkowski, M., Chiavassa, A., Freytag, B., et al. 2016, A&A, 587, A12
  * Wong et al. (2016) Wong, K. T., Kamiński, T., Menten, K. M., & Wyrowski, F. 2016, A&A, 590, A127


# An inverse problem of determining the orders of systems of fractional
pseudo-differential equations

Ravshan Ashurov  and  Sabir Umarov  1Institute of Mathematics of Academy of
Sciences of Republic of Uzbekistan<EMAIL_ADDRESS>2University of New
Haven, Department of Mathematics, 300 Boston Post Road,
West Haven, CT 06516, USA<EMAIL_ADDRESS>
###### Abstract.

As it is known various dynamical processes can be modeled through the systems
of time-fractional order pseudo-differential equations. In the modeling
process one frequently faces with determining adequate orders of time-
fractional derivatives in the sense of Riemann-Liouville or Caputo. This
problem is qualified as an inverse problem. The right (vector) order can be
found utilizing the available data. In this paper we consider this inverse
problem for linear systems of fractional order pseudo-differential equations.
We prove that the Fourier transform of the vector-solution
$\widehat{U}(t,\xi)$ evaluated at a fixed time instance, which becomes
possible due to the available data, recovers uniquely the unknown vector-order
of the system of governing pseudo-differential equations.

Keywords: system of differential equations, fractional order differential
equation, pseudo-differential operator, matrix symbol, inverse problem,
determination of the fractional derivative’s order

##  1\. Introduction

In modern science and engineering researchers frequently use fractional order
differential equations for modeling of dynamics of various complex stochastic
processes arising in different fields; see, for example, [12, 22, 23, 32] in
physics, [20, 28] in finance, [6] in hydrology, [21] in cell biology, among
others. In the last few decades several books, devoted to fractional order
differential and pseudo-differential equations and their various applications,
have been published (see e.g. [27, 24, 16, 11, 29, 31]).

In fractional order modeling, in contrast to integer order equations, orders
of fractional order governing equations are often unknown, and requires to
utilize available data to measure. Therefore, one of the key questions arising
in the process of modeling is a proper determination of a fractional order of
the governing equation. The problem of determining a correct order of an
equation (or orders of equations if the model uses more than one governing
equation) is classified as an inverse problem. Inverse problems naturally
require additional conditions (or information) for a solution. For
subdiffusion equations, in which the order is between zero and one, the
inverse problem of determination of the order has been considered by a number
of authors; see [1, 2, 3, 7, 14, 17, 18, 19] and references therein. For the
survey paper, we refer the reader to [17] by Li et. al. Note that in all the
refereed works the subdiffusion equation was considered in a bounded domain
$\Omega\subset\mathbb{R}^{N}$. In addition, it should be noted that in
publications [7, 14, 18, 19] the following relation was taken as an additional
condition

(1.1) $u(x_{0},t)=h(t),\,\,0<t<T,$

at a monitoring point $x_{0}\in\overline{\Omega}$. But this condition, as a
rule (an exception is the work [14] by J. Janno, where both the uniqueness and
existence are proved), can ensure only the uniqueness of the solution of the
inverse problem [7, 18, 19]. In paper [3] authors obtained the existence and
uniqueness result, considering as an additional condition, the value of the
projection of the solution onto the first eigenfunction of the elliptic part
of the subdiffusion equation. Note that the technique used in [3] is
applicable only when the first eigenvalue is zero. In more general case the
uniqueness and existence of a solution of the inverse problem of determination
of an unknown order of the fractional derivative in the subdiffusion equation
was proved in the recent work [1]. In this case, the additional condition is
$||u(x,t_{0})||^{2}=d_{0}$, and the boundary condition is not necessarily
homogeneous. In paper [2] authors studied the inverse problem for the
simultaneous determination of the order of the Riemann-Liouville time
fractional derivative and the source function in the subdiffusion equations.

The purpose of this work is to investigate the inverse problem of determining
the vector-order of the time-fractional derivatives of systems of pseudo-
differential equations. We note that systems (linear and non-linear) of
fractional ordinary equations and partial differential equations have rich
applications and are used in modeling of various processes arising in modern
science and engineering. For example, they are used in modeling of processes
in biosystems [8, 26, 10], ecology [15, 25], epidemiology [33, 13], etc.

In this paper we consider the following system of linear homogenous time-
fractional order pseudo-differential equations

(1.2)
$\begin{cases}\mathcal{D}^{\beta_{1}}u_{1}(t,x)=A_{1,1}(D)u_{1}(t,x)+\dots
A_{1,m}(D)u_{m}(t,x),\\\
\mathcal{D}^{\beta_{2}}u_{2}(t,x)=A_{2,1}(D)u_{1}(t,x)+\dots
A_{2,m}(D)u_{m}(t,x),\\\ {\cdots}\\\
\mathcal{D}^{\beta_{m}}u_{m}(t,x)=A_{m,1}(D)u_{1}(t,x)+\dots
A_{m,m}(D)u_{m}(t,x).\end{cases}$

where $\mathcal{B}=\langle\beta_{1},\dots,\beta_{m}\rangle,$ $0<\beta_{j}\leq
1,\ j=1,\dots,m,$ is an unknown vector-order to be determined, the operator
$\mathcal{D}$ on the left hand side expresses either the Riemann-Liouville
derivative $D_{+}$ or the Caputo derivative $D_{\ast}$, and $A_{j,k}(D)$ are
pseudo-differential operators with (possibly singular) symbols depending only
on dual variables (for simplicity) and described later. The initial conditions
depend on the form of fractional derivatives.

As it follows from our main result, a predetermined value of the Fourier
transform
$\hat{U}(t,\xi)=\langle\hat{u}_{1}(t,\xi),\dots,\hat{u}_{m}(t,\xi)\rangle$ of
the solution $U(t,x)=\langle u_{1}(t,x),\dots,u_{m}(t,x)\rangle$ of the
initial value problem for system (1.2) at an appropriate fixed point
$\xi_{0}\in\mathbb{R}^{m}$ satisfying some condition (see Eq. (2.11) ), that
is

(1.3) $\hat{u}_{j}(t_{0},\xi^{0})=d_{j},\quad j=1,2,...,m,$

where $t_{0}\geq 1$ is an observation time, uniquely recovers the vector-order
$\mathcal{B}=\langle\beta_{1},\dots,\beta_{m}\rangle$ of the fractional
derivatives.

In the particular case, the determining of a scalar order for one equation was
considered in [3] and the forward problem for systems of pseudo-differential
equations in [30]. From this point of view the current paper is a logical
continuation of these two papers.

We note that in [3] the additional condition is represented in the form

$\int\limits_{\Omega}u(t_{0},x)v_{1}(x)dx=d\neq 0,$

that is, in the form of a projection of the solution $u(t_{0},x)$ onto the
first eigenfunction $v_{1}(x)$ of the elliptic part of the equation considered
in an arbitrary bounded domain $\Omega\subset\mathbb{R}^{n}$. Condition (1.3)
can be considered as the projection of the solution $u(t_{0},x)$ onto ”the
eigenfunction” $e^{-ix\xi_{0}}$:

$\int\limits_{\mathbb{R}^{n}}u(t_{0},x)e^{-ix\xi_{0}}dx=d.$

To our best knowledge, the inverse problem for the system of equations (1.2)
with the vector-order fractional derivatives under the additional condition
(1.3) is considered for the first time.

##  2\. Main results

###  2.1. Notations.

We follow the notions and notations introduced in [30]. For the reader’s
convenience, below we introduce the main notations used in the current paper;
for details see [29, 30]. Let $G\subseteq\mathbb{R}^{n}$ be an open set and
$p\geq 1.$ The space ${\mathbf{\Psi}}_{G,p}(\mathbb{R}^{n})$ comprises of
functions $\psi\in L_{p}(\mathbb{R}^{n})$, such that $\text{supp}\
\hat{\psi}\Subset G$, i.e. the Fourier transform

$\hat{\psi}(\xi)=\int\limits_{\mathbb{R}^{n}}f(x)e^{-ix\xi}dx$

of $\psi$ has a compact support in $G$. This is a topological-vector space
with respect to the following convergence: a sequence $\psi_{n}\to\psi$ if
$\text{supp}\,\hat{\psi}_{n}\Subset G,$ and $\psi_{n}\to\psi$ in
$L_{p}(\mathbb{R}^{n}).$ For relations of the spaces
${\mathbf{\Psi}}_{G,p}(\mathbb{R}^{n})$ to Sobolev spaces and Schwartz
distributions see [29].

Let $A(\xi)$ be a continuous function in $G$. Outside of $G$ or on its
boundary $A(\xi)$ may have singularities of arbitrary type. For a function
$\varphi\in{\mathbf{\Psi}}_{G,p}(\mathbb{R}^{n})$ the pseudo-differential
operator $A(D)$ corresponding to the symbol $A(\xi)$ is defined by the formula

(2.1)
$A(D)\varphi(x)=\frac{1}{(2\pi)^{n}}\int_{G}A(\xi)\hat{\varphi}(\xi)e^{ix\xi}d\xi\,\quad
x\in\mathbb{R}^{n}.$

For the systematic presentation of the theory of pseudo-differential operators
being considered in this paper we refer the reader to [29].

Let

$\mathbb{A}(D)=\begin{bmatrix}A_{1,1}(D)&\dots&A_{1,m}(D)\\\
\dots&\dots&\dots\\\ A_{m,1}(D)&\dots&A_{m,m}(D)\end{bmatrix}$

be the matrix pseudo-differential operator with constant (that is not
depending on the variable $x$) matrix-symbol

(2.2) $\mathcal{A}(\xi)=\begin{bmatrix}A_{1,1}(\xi)&\dots&A_{1,m}(\xi)\\\
\dots&\dots&\dots\\\ A_{m,1}(\xi)&\dots&A_{m,m}(\xi)\end{bmatrix},$

defined and continuous in $G$ in the sense of the matrix norm.

With the matrix form of the pseudo-differential operator we can represent
system (1.2) in the vector form:

(2.3) $\mathcal{D}^{\mathcal{B}}{U}(t,x)=\mathbb{A}(D){U}(t,x),$

where
$\mathcal{D}^{\mathcal{B}}{U}(t,x)=\langle\mathcal{D}^{{\beta_{1}}}{u}_{1}(t,x),\dots,\mathcal{D}^{{\beta_{m}}}{u}_{m}(t,x)\rangle.$

Below we will use two main forms of fractional derivatives, namely the
Riemann-Liouville form and the Caputo form. Let $k$ be a natural number and
$k-1\leq\beta<k.$ Then the fractional derivative of order $\beta$ of a
measurable function $f$ in the sense of Riemann–Liouville is defined as

$D^{\beta}_{+}f(t)=\frac{1}{\Gamma(k-\beta)}\frac{d^{k}}{dt^{k}}\int_{0}^{t}\frac{f(\tau)d\tau}{(t-\tau)^{\beta+1-k}},$

provided the expression on the right exists. Here $\Gamma(t)$ is Euler’s gamma
function. If we replace differentiation and fractional integration in this
definition, then we get the definition of a regularized derivative, that is,
the definition of a fractional derivative in the sense of Caputo:

$D^{\beta}_{\ast}f(t)=\frac{1}{\Gamma(k-\beta)}\int_{0}^{t}\frac{f^{(k)}(\tau)d\tau}{(t-\tau)^{\beta+1-k}},$

provided the integral on the right exists.

We assume that the matrix-symbol is symmetric, $A_{k,j}(\xi)=A_{j,k}(\xi)$ for
all $k,j=1,\dots,m,$ and $\xi\in G,$ and diagonalizable. Namely, there exists
an invertible $(m\times m)$-matrix-function $M(\xi),$ such that

(2.4) $\mathcal{A}(\xi)=M^{-1}(\xi)\Lambda(\xi)M(\xi),\quad\xi\in G,$

with a diagonal matrix

(2.5) ${\Lambda}(\xi)=\begin{bmatrix}\lambda_{1}(\xi)&\dots&0\\\
\dots&\dots&\dots\\\ 0&\dots&\lambda_{m}(\xi)\end{bmatrix}.$

We denote entries of matrices $M(\xi)$ and $M^{-1}(\xi)$ by $\mu_{j,k}(\xi),\
j,k=1,\dots,m,$ and $\nu_{j,k}(\xi)$, $j,k=1,\dots,m,$ respectively.

Since initial conditions depend on the form of the fractional derivative on
the left hand side of equation (2.3), we will consider the cases with the
Caputo and Riemann-Liouville derivatives separately. We first formulate our
main result in the case of Caputo fractional derivative. The case of Rieman-
Liouville fractional derivative can be treated similarly.

###  2.2. Forward problem

Let $\mathcal{B}$ be a known vector-order with $0<\beta_{j}\leq 1$,
$j=1,\dots,m$. Consider the following Cauchy problem

(2.6) $\displaystyle D_{\ast}^{\mathcal{B}}{U}(t,x)$
$\displaystyle=\mathbb{A}(D){U}(t,x),\quad t>0,\ x\in\mathbb{R}^{n},$ (2.7)
$\displaystyle U(0,x)$ $\displaystyle=\varPhi(x),\quad x\in\mathbb{R}^{n},$

where
$\varPhi(x)=\langle\varphi_{1}(x),\dots,\varphi_{m}(x)\rangle\in{\mathbf{\Psi}}_{G,p}(\mathbb{R}^{n})$
and the fractional derivatives on the left are in the sense of Caputo.

We call the Cauchy problem (2.6)-(2.7) _the forward problem._

A representation formula for the solution of the forward problem was obtained
in [30] and it has the form

$u_{j}(t,x)=\frac{1}{(2\pi)^{n}}\int\limits_{\mathbb{R}^{n}}\sum\limits_{k=1}^{m}s_{j,k}(t,\xi)\hat{\varphi}_{k}(\xi)e^{ix\xi}d\xi,\quad
j=1,\dots,m,$

where

$s_{j,k}(t,\xi)=\sum\limits_{l=1}^{m}\mu_{j,l}(\xi)\nu_{l,k}(\xi)E_{\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}}).$

Here we denoted by $E_{\beta_{j}}(z),j=1,\dots,m,$ the Mittag-Leffler
functions of indices $\beta_{1},\dots,\beta_{m},$ respectively.

We rewrite function $u_{j}$ as

$\displaystyle u_{j}(t,x)$
$\displaystyle=\frac{1}{(2\pi)^{n}}\int\limits_{\mathbb{R}^{n}}\sum\limits_{k=1}^{m}\sum\limits_{l=1}^{m}\mu_{j,l}(\xi)\nu_{l,k}(\xi)E_{\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}})\hat{\varphi}_{k}(\xi)e^{ix\xi}d\xi$
$\displaystyle=\frac{1}{(2\pi)^{n}}\int\limits_{\mathbb{R}^{n}}\sum\limits_{l=1}^{m}E_{\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}})\bigg{[}\mu_{j,l}(\xi)\sum\limits_{k=1}^{m}\nu_{l,k}(\xi)\hat{\varphi}_{k}(\xi)\bigg{]}e^{ix\xi}d\xi$
$\displaystyle=\frac{1}{(2\pi)^{n}}\int\limits_{\mathbb{R}^{n}}\sum\limits_{l=1}^{m}E_{\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}})K_{j,l}\big{(}\xi,\widehat{\Phi}(\xi)\big{)}e^{ix\xi}d\xi,$

where

$K_{j,l}\big{(}\xi,\widehat{\Phi}(\xi)\big{)}=\mu_{j,l}(\xi)\sum\limits_{k=1}^{m}\nu_{l,k}(\xi)\hat{\varphi}_{k}(\xi),$

and

$\widehat{\Phi}(\xi)=\langle\hat{\varphi}_{1}(\xi),\hat{\varphi}_{2}(\xi),\cdot\cdot\cdot,\hat{\varphi}_{m}(\xi)\rangle.$

For the Fourier transform of the solution we have

(2.8)
$\hat{u}_{j}(t,\xi)=\sum\limits_{l=1}^{m}E_{\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}})K_{j,l}\big{(}\xi,\widehat{\Phi}(\xi)\big{)}.$

Note that under the above conditions on the matrix-symbol $A_{j,k}(\xi)$ and
on the function $\Phi(x)$, this Fourier transform exists at each point
$\xi\in\mathbb{R}^{n}$.

###  2.3. Inverse problem.

Now let the parameter $\mathcal{B}$ be an unknown vector-order of the time
derivative with $\beta_{0}\leq\beta_{j}<1$, $j=1,\dots,m$,
$\beta_{0}\in(0,1)$. The main purpose of this paper is to investigate the
inverse problem of identifying of these parameters $\beta_{j}$. Since there
are $m$ unknown parameters, we need to set $m$ conditions. We pass on to the
determining of these additional conditions.

In what follows, we will assume that

(2.9) $|\arg\lambda_{j}(\xi)|>\frac{\pi}{2},\quad\xi\in G,\quad j=1,\dots,m.$

Let $\xi^{0}=(\xi^{0}_{1},\xi^{0}_{2},...,\xi^{0}_{m})\in
G\subset\mathbb{R}^{n}$ be a vector such that the determinant of the matrix

(2.10)
$\mathcal{K}(\xi^{0})\equiv\\{K_{j,l}\big{(}\xi^{0},\widehat{\Phi}(\xi^{0})\big{)}\\},\
j,l=1,\dots,m,$

satisfies the condition

(2.11) $\big{|}K_{j,l}\big{(}\xi^{0},\widehat{\Phi}(\xi^{0})\big{)}\big{|}\neq
0.$

To find the unknown parameters $\beta_{l},l=1,\dots,m$, we consider the
following additional conditions

(2.12)
$f_{j}(\mathcal{B},t_{0},\xi^{0})\equiv\hat{u}_{j}(t_{0},\xi^{0})=d_{j},\quad
j=1,\dots,m,$

where $d_{j}$ are given numbers and $t_{0}$ is defined later (see Lemmas 2.2
and 2.4). We call problem (2.6)–(2.7) together with the additional condition
(2.12) _the inverse problem._

It follows from (2.8) and (2.12) that for all $j=1,2,...,m$

(2.13)
$\sum\limits_{l=1}^{m}E_{\beta_{l}}(\lambda_{l}(\xi^{0})t_{0}^{\beta_{l}})K_{j,l}\big{(}\xi^{0},\widehat{\Phi}(\xi^{0})\big{)}=d_{j}.$

These are in fact the system of equations to define the orders
$\beta_{l},l=1,2,...,m$.

Due to condition (2.11) one can solve system (2.13) with respect to the
Mittag-Leffler functions $E_{\beta_{l}}$, i.e.

(2.14) $E_{\beta_{l}}(\lambda_{l}(\xi^{0})t_{0}^{\beta_{l}})=b_{l},\quad
l=1,2,...,m,$

where $b_{l},\,l=1,2,...,m,$ are components of the vector
$\mathcal{K}^{-1}(\xi^{0})\textbf{d},$ with $\textbf{d}=\langle
d_{1},\dots,d_{m}\rangle$ and $\mathcal{K}^{-1}(\xi^{0})$ is the inverse
matrix to $\mathcal{K}(\xi^{0}),$ defined in (2.10). Thus to define each
unknown parameter $\beta_{l}$ we obtained a separate equation (2.14).

Let $R_{C,l}$ be the range of values of the function
$e_{1,\lambda_{l}}(\beta_{l})\equiv
E_{\beta_{l}}(\lambda_{l}(\xi^{0})t_{0}^{\beta_{l}})$ when $\beta_{l}$ runs
over the half-interval $[\beta_{0},1)$, i.e.

$e_{1,\lambda_{l}}:[\beta_{0},1)\rightarrow R_{C,l}\subset C,$

where $C$ is a complex plane, and the index $C$ emphasizes that we are
considering the case of the Caputo derivatives.

Obviously, for equations (2.14), in order to have solutions with respect to
$\beta_{l},$ the right-hand sides of these equations must lie within the
values of the functions on the left-hand sides of these equations, i.e.

(2.15) $b_{l}\in R_{C,l},\quad l=1,2,...,m.$

On the other hand, by virtue of Rolle’s theorem, the strict monotonicity of
either the real part or the imaginary part of the function
$E_{\beta_{l}}(\lambda_{l}(\xi^{0})t_{0}^{\beta_{l}})$ in the variable
$\beta_{l}$ is sufficient for the uniqueness of the solution to equation
(2.14).

Let us introduce the following notation

$R_{C}(\beta_{l})=\Re(E_{\beta_{l}}(\lambda_{l}(\xi^{0})t_{0}^{\beta_{l}})),$

where $\Re(z)$ is the real part of $z$ and the index $C$ again emphasizes that
we are considering the case of the Caputo derivatives. The necessity of
condition (2.9) is that its fulfillment guarantees the strict monotonicity of
the function $R_{C}(\beta_{l})$ it the variable $\beta_{l}$ for each fixed
$l.$

###  2.4. Main results

The main results of this paper are stated in Theorems 2.1 and 2.5.

######  Theorem 2.1.

Let $\xi^{0}$ satisfy condition (2.11) and $t_{0}>T_{0},$ where $T_{0}$ is
identified in Lemma 2.2. Let the numbers $d_{l}$ on the right hand side of
equation (2.12) be such that the components $b_{l},\ l=1,\dots,m,$ of the
vector $\mathcal{K}^{-1}(\xi^{0})\textbf{d}$ satisfy the conditions (2.15).
Then for each $l$ there exists the unique number
$\beta_{l}^{\ast}\in[\beta_{0},1]$ such that the Fourier transform of the
solution $u_{j}(t,x),j=1,\dots,m,$ of the forward problem with
$\beta_{j}=\beta_{j}^{\ast},j=1,\dots,m,$ satisfies equation (2.12).

The proof of this theorem follows from the existence and uniqueness theorem
for the forward problem proved in [30] (see Theorem 3.1) and Lemma 2.2 below.
Therefore, in order to prove the theorem we need to prove only this lemma. The
proof of Lemma 2.2 is given in Section 3.

######  Lemma 2.2.

Given $\beta_{0}$ in the interval $0<\beta_{0}<1$ and $\xi_{0}$ satisfying
condition (2.11), there exists a number $T_{0}=T_{0}(\xi^{0},\beta_{0})>1$,
such that for all $t_{0}\geq T_{0}$ the function $R_{C}(\beta_{l})$ is
positive and strictly monotonically decreasing with respect to
$\beta_{l}\in[\beta_{0},1]$ and

(2.16) $R_{C}(1)\leq R_{C}(\beta_{l})\leq R_{C}(\beta_{0}),\,\,l=1,\dots,m.$

######  Remark 2.3.

Theorem 2.1 defines the vector-order
$\mathcal{B}^{\ast}=(\beta_{1}^{\ast},\beta_{2}^{\ast},...,\beta_{m}^{\ast})$
uniquely from conditions (2.12). Hence, if we define
$f_{j}(\mathcal{B},\cdot,\cdot)$ at another time instant $t_{1}$ and point
$\xi^{1}$ and get a new $\mathcal{B^{\ast\ast}}$, i.e.
$f_{j}(\mathcal{B}^{\ast\ast},t_{1},\xi^{1})=d_{j}^{1}$, then from the
equality $f_{j}(\mathcal{B}^{\ast\ast},t_{0},\xi^{0})=d_{j}$, by virtue of the
theorem, we obtain $\mathcal{B}^{\ast\ast}=\mathcal{B}^{\ast}$.

Now consider the following initial-value problem

(2.17) $\displaystyle D_{+}^{\mathcal{B}}{U}(t,x)$
$\displaystyle=\mathbb{A}(D){U}(t,x),\quad t>0,\ x\in\mathbb{R}^{n},$ (2.18)
$\displaystyle J^{1-\mathcal{B}}U(0,x)$ $\displaystyle=\varPhi(x),\quad
x\in\mathbb{R}^{n},$

where
$\varPhi(x)=\langle\varphi_{1}(x),\dots,\varphi_{m}(x)\rangle\in{\mathbf{\Psi}}_{G,p}(\mathbb{R}^{n})$
and the fractional derivatives on the left hand side of equation (2.17) are in
the sense of Riemann-Liouville.

We call the Cauchy problem (2.17)-(2.18) _the second forward problem._

A representation formula for the solution of the second forward problem was
also obtained in [30] and it has the form

$u_{j}(t,x)=\frac{1}{(2\pi)^{n}}\int\limits_{\mathbb{R}^{n}}\sum\limits_{k=1}^{m}s^{+}_{j,k}(t,\xi)\hat{\varphi}_{k}(\xi)e^{ix\xi}d\xi,\quad
j=1,\dots,m,$

where

$s^{+}_{j,k}(t,\xi)=\sum\limits_{l=1}^{m}\mu_{j,l}(\xi)\nu_{l,k}(\xi)t^{\beta_{l}-1}E_{\beta_{l},\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}}).$

Here we denoted by $E_{\beta_{j},\beta_{j}}(z),j=1,\dots,m,$ the two-
parametric Mittag-Leffler functions.

For the Fourier transform of the solution we have

(2.19)
$\hat{u}_{j}(t,\xi)=\sum\limits_{l=1}^{m}t^{\beta_{l}-1}E_{\beta_{l},\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}})K_{j,l}\big{(}\xi,\widehat{\Phi}(\xi)\big{)}.$

Suppose that condition (2.9) is fulfilled and choose $\xi^{0}\in G$ so that
inequality (2.11) holds.

We call problem (2.17)-(2.18) together with the additional condition (2.12)
_the second inverse problem_.

Note that additional condition (2.12) is in fact the equation to determine the
unknown parameters $\beta_{l}$. Performing similar calculations as above, by
virtue of condition (2.11), we rewrite (2.12) as

(2.20)
$t^{\beta_{l}-1}E_{\beta_{l},\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}})=b_{l},\quad
l=1,2,...,m,$

where $b_{l},\,l=1,2,...,m,$ are the same numbers as above.

Let $R_{RL,l}$ be the range of values of the left-hand side of these equations
when $\beta_{l}$ runs over the half-interval $[\beta_{0},1)$. Here the index
$RL$ emphasizes that we are considering the case of the Riemann-Liouville
derivatives.

Again, as in case of equations (2.14), a necessary condition for the existence
of solutions to equations (2.20) is the inclusion

(2.21) $b_{l}\in R_{RL,l},\quad l=1,2,...,m.$

On the other hand, by virtue of Rolle’s theorem, the strict monotonicity of
the function
$\Re(t^{\beta_{l}-1}E_{\beta_{l},\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}}))$
in the variable $\beta_{l}$ is sufficient for the uniqueness of the solution
to equation (2.20).

However, if $|\Re(\lambda_{l}(\xi^{0}))|=|\Im(\lambda_{l}(\xi^{0}))|$ (note,
under the condition (2.9) one has $\Re(\lambda_{l}(\xi^{0}))<0$), then the
principal part of
$\Re(t^{\beta_{l}-1}E_{\beta_{l},\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}}))$
vanishes, and in this case it is necessary to go to its next term in the
asymptotic. Therefore, to simplify the presentation, we further assume that

(2.22) $|\Re(\lambda_{l}(\xi^{0}))|\neq|\Im(\lambda_{l}(\xi^{0}))|$

Let us introduce the following notation

$R_{RL}(\beta_{l})=sign(|\Re(\lambda_{l}(\xi^{0}))|-|\Im(\lambda_{l}(\xi^{0}))|)\Re(t^{\beta_{l}-1}E_{\beta_{l},\beta_{l}}(\lambda_{l}(\xi)t^{\beta_{l}})).$

Here the index $RL$ again emphasizes that we are considering the case of the
Riemann-Liouville derivatives.

######  Lemma 2.4.

Given $\beta_{0}$ in the interval $0<\beta_{0}<1$ and $\xi_{0}$ satisfying
condition (2.11), there exists a number $T_{1}=T_{1}(\xi^{0},\beta_{0})>1$,
such that for all $t_{0}\geq T_{1}$ the function $R_{RL}(\beta_{l})$ is
positive and strictly monotonically decreasing with respect to
$\beta_{l}\in[\beta_{0},1]$ and

(2.23) $R_{RL}(1)\leq R_{RL}(\beta_{l})\leq
R_{RL}(\beta_{0}),\,\,l=1,\dots,m.$

This lemma, similar to the Caputo derivative case, immediately implies the
following main result of this paper in the case of the Riemann-Liouville
derivatives. The existence and uniqueness theorem of the corresponding forward
problem is proved in [30] (see Theorem 3.4). The proof of Lemma 2.4 is
presented in the next section.

######  Theorem 2.5.

Let $\xi^{0}$ satisfy condition (2.11) and $t_{0}>T_{1},$ where $T_{1}$ is
identified in Lemma 2.4. Let the numbers $d_{l}$ from condition (2.12) be such
that the corresponding numbers $b_{l}$ satisfy the conditions (2.21). Then for
each $l$ there exists the unique number $\beta_{l}^{\ast}\in[\beta_{0},1]$
such that the Fourier transform of the solution $u_{j}(t,x)$ of the second
forward problem with $\beta_{j}=\beta_{j}^{\ast}$ satisfies the equation
(2.12).

Similar to the Caputo derivative case, Theorem 2.5 defines the vector-order
$\mathcal{B}^{\ast}=(\beta_{1}^{\ast},\beta_{2}^{\ast},...,\beta_{m}^{\ast})$
uniquely from conditions (2.12); see Remark 2.3.

##  3\. Proofs of Lemmata 2.2 and 2.4

Let us denote by $\delta(1;\theta)$ a contour oriented by non-decreasing
$\arg\zeta$ consisting of the following parts: the ray $\arg\zeta=-\theta$
with $|\zeta|\geq 1$, the arc $-\theta\leq\arg\zeta\leq\theta$, $|\zeta|=1$,
and the ray $\arg\zeta=\theta$, $|\zeta|\geq 1$. If $0<\theta<\pi$, then the
contour $\delta(1;\theta)$ divides the complex $\zeta$-plane into two
unbounded parts, namely $G^{(-)}(1;\theta)$ to the left of $\delta(1;\theta)$
by orientation, and $G^{(+)}(1;\theta)$ to the right of it. The contour
$\delta(1;\theta)$ is called the Hankel path.

In what follows, we fix $l$ out of $1,\dots,m$ and denote
$\lambda=\lambda_{l}(\xi^{0})$ and $\rho=\beta_{l}$. Let
$\lambda=-\lambda_{1}+i\lambda_{2}$ and by virtue of condition (2.9) one has
$\lambda_{1}>0$. Further let $\theta=(\frac{\pi}{2}+\varepsilon)\rho$,
$\alpha=(\frac{\pi}{2}+2\varepsilon)\rho$, $\rho\in[\beta_{0},1)$ and
$\varepsilon>0$ be such that
$\varepsilon\equiv\varepsilon(\xi^{0})<\frac{1}{2}\min\\{|\arg\lambda(\xi^{0})|-\pi/2,\pi/2\\}$.

Then

$\frac{\pi}{2}\rho<\theta<\alpha<\pi\rho,\quad\alpha<|\arg\lambda|,$

and therefore $\lambda t_{0}^{\rho}\in G^{(-)}(1;\theta)$.

###  3.1. Proof of Lemma 2.2

First we prove Lemma 2.2. By the definition of contour $\delta(1;\theta)$, we
have (see [9], formula (2.29), p. 135)

(3.1) $e_{1,\lambda}(\rho)\equiv E_{\rho}(\lambda
t_{0}^{\rho})=-\frac{1}{\lambda t_{0}^{\rho}\Gamma(1-\rho)}+\frac{1}{2\pi
i\rho\lambda
t_{0}^{\rho}}\int\limits_{\delta(1;\theta)}\frac{e^{\zeta^{1/\rho}}\zeta}{\zeta+\lambda
t_{0}^{\rho}}d\zeta=f_{1}(\rho)+f_{2}(\rho).$

To prove the lemma, we need to determine the sign of the real part of the
derivative $\frac{d}{d\rho}e_{1,\lambda}(\rho)$. It is not hard to estimate
the derivative $f^{\prime}_{1}(\rho)$. Indeed, let $\Psi(\rho)$ be the
logarithmic derivative of the gamma function $\Gamma(\rho)$ (for the
definition and properties of $\Psi$ see [5]). Then
$\Gamma^{\prime}(\rho)=\Gamma(\rho)\Psi(\rho)$, and therefore,

$f_{1}^{\prime}(\rho)=\frac{\ln t_{0}-\Psi(1-\rho)}{\lambda
t_{0}^{\rho}\Gamma(1-\rho)}.$

Since

$\frac{1}{\Gamma(1-\rho)}=\frac{1-\rho}{\Gamma(2-\rho)},\quad\Psi(1-\rho)=\Psi(2-\rho)-\frac{1}{1-\rho},$

the function $f_{1}^{\prime}(\rho)$ can be represented as follows

$f_{1}^{\prime}(\rho)=\frac{1}{\lambda t_{0}^{\rho}}\frac{(1-\rho)[\ln
t_{0}-\Psi(2-\rho)]+1}{\Gamma(2-\rho)}.$

If $\gamma\approx 0,57722$ is the Euler-Mascheroni constant, then
$-\gamma<\Psi(2-\rho)<1-\gamma$. By virtue of this estimate we may write

(3.2)
$-\Re(f_{1}^{\prime}(\rho))\geq\frac{\lambda_{1}}{|\lambda|^{2}}\frac{(1-\rho)[\ln
t_{0}-(1-\gamma)]+1}{\Gamma(2-\rho)t_{0}^{\rho}}\geq\frac{\lambda_{1}}{|\lambda|^{2}t_{0}^{\rho}},$

provided

(3.3) $\ln t_{0}>1-\gamma\quad\text{or}\quad t_{0}>T_{0}=e^{1-\gamma}>1.$

To estimate the derivative $f^{\prime}_{2}(\rho)$, we denote the integrand in
(3.1) by $F(\zeta,\rho)$:

$F(\zeta,\rho)=\frac{1}{2\pi i\rho\lambda
t_{0}^{\rho}}\cdot\frac{e^{\zeta^{1/\rho}}\zeta}{\zeta+\lambda t_{0}^{\rho}}.$

Note, that the domain of integration $\delta(1;\theta)$ also depends on
$\rho$. To take this circumstance into account when differentiating the
function $f^{\prime}_{2}(\rho)$, we rewrite the integral (3.1) in the form:

$f_{2}(\rho)=f_{2+}(\rho)+f_{2-}(\rho)+f_{21}(\rho),$

where

$f_{2\pm}(\rho)=e^{\pm i\theta}\int\limits_{1}^{\infty}F(s\,e^{\pm
i\theta},\rho)\,ds,$

$f_{21}(\rho)=i\int\limits_{-\theta}^{\theta}F(e^{iy},\rho)\,e^{iy}dy=i\theta\int\limits_{-1}^{1}F(e^{i\theta
s},\rho)\,e^{i\theta s}ds.$

Let us consider the function $f_{2+}(\rho)$. Since
$\theta=(\frac{\pi}{2}+\varepsilon)\rho$ and $\zeta=s\,e^{i\theta}$, then

$e^{\zeta^{1/\rho}}=e^{-s^{\frac{1}{\rho}}(\varepsilon_{1}-i\varepsilon_{2})},\,\,\cos(\frac{\pi}{2}+\varepsilon)=-\varepsilon_{1}<0.\,\,\sin(\frac{\pi}{2}+\varepsilon)=\varepsilon_{2}>0.$

The derivative of the function $f_{2+}(\rho)$ has the form

$f_{2+}^{\prime}(\rho)=\frac{1}{2\pi i\rho\lambda
t_{0}^{\rho}}\int\limits_{1}^{\infty}\frac{e^{-s^{\frac{1}{\rho}}(\varepsilon_{1}-i\varepsilon_{2})}\,s\,e^{2ia\rho}\mathcal{M}(s)}{s\,e^{ia\rho}+\lambda
t_{0}^{\rho}}ds,$

where $a=\frac{\pi}{2}+\varepsilon,$ and

$\mathcal{M}(s)=-\frac{\varepsilon_{1}-i\varepsilon_{2}}{\rho^{2}}s^{1/\rho}\ln
s+2ia-\frac{1}{\rho}-\ln t_{0}-\frac{ias\,e^{ia\rho}+\lambda t_{0}^{\rho}\ln
t_{0}}{s\,e^{ia\rho}+\lambda t_{0}^{\rho}}.$

It is not hard to verify, that

$|s\,e^{ia\rho}+\lambda
t_{0}^{\rho}|\geq|\lambda|t_{0}^{\rho}\sin(\alpha-\theta)\geq\frac{2}{\pi}|\lambda|t_{0}^{\rho}\varepsilon.$

Therefor we arrive at

$|f_{2+}^{\prime}(\rho)|\leq\frac{C}{\rho(\varepsilon|\lambda|t_{0}^{\rho})^{2}}\int\limits_{1}^{\infty}e^{-\varepsilon_{1}\,s^{1/\rho}}s\,\left[\frac{1}{\rho^{2}}s^{1/\rho}\ln
s+\ln t_{0}\right]ds,$

or

$|f_{2+}^{\prime}(\rho)|\leq\frac{C}{(|\lambda|t_{0}^{\rho})^{2}}\left[\frac{1}{\rho}+\ln
t_{0}\right],$

where the constant $C$ depends only on $\varepsilon$ (and therefore only on
$\xi^{0}$).

The function $f^{\prime}_{2-}(\rho)$ has exactly the same estimate.

Now consider the function $f_{21}(\rho)$. It is not hard to verify that

$f^{\prime}_{21}(\rho)=\frac{a}{2\pi\lambda
t_{0}^{\rho}}\int\limits_{-1}^{1}\frac{e^{e^{ias}}e^{2ia\rho s}\big{[}2ias-\ln
t_{0}-\frac{iase^{ia\rho s}+\lambda t_{0}^{\rho}\ln t_{0}}{e^{ia\rho
s}+\lambda t_{0}^{\rho}}\big{]}}{e^{ia\rho s}+\lambda t_{0}^{\rho}}ds.$

Therefore,

$|f^{\prime}_{21}(\rho)|\leq C\frac{\ln t_{0}}{(|\lambda|t_{0}^{\rho})^{2}}.$

Now we show that the real part of the derivative
$\frac{d}{d\rho}e_{1,\lambda}(\rho)$ is negative. Taking into account estimate
(3.2) and the estimates for $f^{\prime}_{2\pm}$ and $f^{\prime}_{21}$, we have

(3.4)
$\Re\bigg{(}\frac{d}{d\rho}e_{1,\lambda}(\rho)\bigg{)}<-\frac{\lambda_{1}}{|\lambda|^{2}t_{0}^{\rho}}+C\frac{1/\rho+\ln
t_{0}}{(|\lambda|t_{0}^{\rho})^{2}}.$

In other words, this derivative is negative if

$t_{0}^{\rho}>C\frac{1/\rho+\ln t_{0}}{\lambda_{1}}$

for all $\rho\in[\beta_{0},1).$ Hence

(3.5) $t_{0}^{\beta_{0}}>C\frac{1/{\beta_{0}}+\ln t_{0}}{\lambda_{1}}.$

Thus, there exists a number $T_{0}=T_{0}(\xi^{0},\beta_{0})>1$ (see (3.3))
such, that for all $t_{0}\geq T_{0}$ we have the estimate

$\Re\bigg{(}\frac{d}{d\rho}e_{1,\lambda}(\rho)\bigg{)}<0\,\,\text{for}\,\,\text{all}\quad\rho\in[\beta_{0},1].$

The positivity of $R_{C}(\beta_{l})$ follows from the explicit form of the
function $f_{1}(\rho)$.

Lemma 2.2, and therefore Theorem 2.1 are completely proved.

###  3.2. Proof of Lemma 2.4

We now turn to the proof of Lemma 2.4. By the definition of contour
$\delta(1;\theta)$, we have for $e_{2,\lambda}(\rho)\equiv
t_{0}^{\rho-1}E_{\rho,\rho}(\lambda t_{0}^{\rho})$ the following equation (see
[9], formula (2.29), p. 135)

(3.6)
$e_{2,\lambda}(\rho)=-\frac{1}{\lambda^{2}t_{0}^{\rho+1}\Gamma(-\rho)}+\frac{1}{2\pi
i\rho\lambda^{2}t_{0}^{\rho+1}}\int\limits_{\delta(1;\theta)}\frac{e^{\zeta^{1/\rho}}\zeta^{\frac{1}{\rho}+1}}{\zeta+\lambda
t_{0}^{\rho}}d\zeta=g_{1}(\rho)+g_{2}(\rho).$

Since the positivity of $R_{RL}(1)$ is obvious, then in order to prove Lemma
2.4 it suffices to show that the derivatives of $R_{RL}(\rho)$ is negative for
all $\rho\in[\beta_{0},1)$.

For the derivative $g^{\prime}_{1}(\rho)$ we have

$g_{1}^{\prime}(\rho)=\frac{\ln
t_{0}-\Psi(-\rho)}{\lambda^{2}t_{0}^{\rho+1}\Gamma(-\rho)}.$

To get rid of the singularity in the denominators, we use the equalities

$\displaystyle\frac{1}{\Gamma(-\rho)}$
$\displaystyle=-\frac{\rho}{\Gamma(1-\rho)}=-\frac{\rho(1-\rho)}{\Gamma(2-\rho)},$
$\displaystyle\Psi(-\rho)$
$\displaystyle=\Psi(1-\rho)+\frac{1}{\rho}=\Psi(2-\rho)+\frac{1}{\rho}-\frac{1}{1-\rho}.$

Then the function $g_{1}^{\prime}(\rho)$ can be represented as follows

(3.7)
$g_{1}^{\prime}(\rho)=\frac{1}{\lambda^{2}t_{0}^{\rho+1}}\frac{\rho(1-\rho)[\Psi(2-\rho)-\ln
t_{0}]+1-2\rho}{\Gamma(2-\rho)}=-\frac{g_{11}(\rho)}{\lambda^{2}t_{0}^{\rho+1}\Gamma(2-\rho)}.$

Since $\Psi(2-\rho)<1-\gamma$, then

$g_{11}(\rho)>\rho(1-\rho)[\ln t_{0}-(1-\gamma)])+2\rho-1.$

For $t_{0}=e^{1-\gamma}e^{2/\rho}$ one has $\rho(1-\rho)[\ln
t_{0}-(1-\gamma)])+2\rho-1=1$. Hence, $g_{11}(\rho)\geq 1$, provided
$t_{0}\geq T_{1},$ where

(3.8) $T_{1}=e^{1-\gamma}e^{2/\beta_{0}}>e^{3-\gamma}>1.$

Thus, by virtue of (3.7), for all such $t_{0}$ we arrive at

(3.9)
$sign(\lambda_{1}^{2}-\lambda_{2}^{2})\Re(g_{1}^{\prime}(\rho))\leq-\frac{|\lambda_{1}^{2}-\lambda_{2}^{2}|}{|\lambda|^{4}t_{0}^{\rho+1}}.$

To estimate the derivative $g^{\prime}_{2}(\rho)$, we denote the integrand in
(3.6) by $G(\zeta,\rho)$:

$G(\zeta,\rho)=\frac{1}{2\pi
i\rho\lambda^{2}t_{0}^{\rho+1}}\cdot\frac{e^{\zeta^{1/\rho}}\zeta^{1/\rho+1}}{\zeta+\lambda
t_{0}^{\rho}},$

and rewrite the integral (3.6) in the form:

$g_{2}(\rho)=g_{2+}(\rho)+g_{2-}(\rho)+g_{21}(\rho),$

where

$g_{2\pm}(\rho)=e^{\pm i\theta}\int\limits_{1}^{\infty}G(s\,e^{\pm
i\theta},\rho)\,ds,$

$g_{21}(\rho)=i\int\limits_{-\theta}^{\theta}F(e^{iy},\rho)\,e^{iy}dy=i\theta\int\limits_{-1}^{1}G(e^{i\theta
s},\rho)\,e^{i\theta s}ds.$

The derivative of the function $g_{2+}(\rho)$ has the form

$g_{2+}^{\prime}(\rho)=I\cdot\int\limits_{1}^{\infty}\frac{e^{-s^{1/\rho}\,(\varepsilon_{1}-i\varepsilon_{2})}s^{\frac{1}{\rho}+1}\,e^{2ia\rho}\mathcal{N}(s)}{s\,e^{ia\rho}+\lambda
t_{0}^{\rho}}ds,$

where $I=e^{ia}(2\pi i\rho\lambda^{2}t_{0}^{\rho+1})^{-1},$
$a=\frac{3\pi}{4},$ and

$\mathcal{N}(s)=-\frac{1}{\rho^{2}}((\varepsilon_{1}+i\varepsilon_{2})s^{\frac{1}{\rho}}+1)\ln
s+2ia-\frac{1}{\rho}-\ln t_{0}-\frac{ias\,e^{ia\rho}+\lambda t_{0}^{\rho}\ln
t_{0}}{s\,e^{ia\rho}+\lambda t_{0}^{\rho}}.$

By virtue of the inequality $|s\,e^{ia\rho}+\lambda
t_{0}^{\rho}|\geq\frac{2}{\pi}|\lambda|t_{0}^{\rho}\varepsilon$ we arrive at

$|g_{2+}^{\prime}(\rho)|\leq\frac{C}{\rho|\lambda|^{3}t_{0}^{2\rho+1}}\int\limits_{1}^{\infty}e^{-\varepsilon_{1}\,s^{1/\rho}}s^{\frac{1}{\rho}+1}\,\big{[}\frac{1}{\rho^{2}}s^{1/\rho}\ln
s+\ln t_{0}\big{]}ds,$

or

$|g_{2+}^{\prime}(\rho)|\leq\frac{C}{|\lambda|^{3}t_{0}^{2\rho+1}}\,\big{[}\frac{1}{\rho}+\ln
t_{0}\big{]},$

where the constant $C$ depends only on $\xi^{0}$.

The function $g^{\prime}_{2-}(\rho)$ has exactly the same estimate.

Now consider the function $g_{21}(\rho)$. For its derivative we have

$g^{\prime}_{21}(\rho)=\frac{a}{2\pi
i\lambda^{2}t_{0}^{\rho+1}}\cdot\int\limits_{-1}^{1}\frac{e^{e^{ias}}\,e^{ias}\,e^{2ia\rho
s}\big{[}2ias-\ln t_{0}-\frac{iase^{ia\rho s}+\lambda t_{0}^{\rho}\ln
t_{0}}{e^{ia\rho s}+\lambda t_{0}^{\rho}}\big{]}}{e^{ia\rho s}+\lambda
t_{0}^{\rho}}ds.$

Therefore,

$|g^{\prime}_{21}(\rho)|\leq C\,\frac{\ln
t_{0}}{|\lambda|^{3}t_{0}^{2\rho+1}}.$

Taking into account estimate (3.9) and the estimates for $g^{\prime}_{2\pm}$
and $g^{\prime}_{21}$, we have

$sign(\lambda_{1}^{2}-\lambda_{2}^{2})\Re(e^{\prime}_{2,\lambda}(\rho))\leq-\frac{|\lambda_{1}^{2}-\lambda_{2}^{2}|}{|\lambda|^{4}t_{0}^{\rho+1}}+C\frac{1/\rho+\ln
t_{0}}{|\lambda|^{3}t_{0}^{2\rho+1}}.$

In other words, the left hand side is negative if

$t_{0}^{\beta_{0}}>\frac{C|\lambda|}{|\lambda_{1}^{2}-\lambda_{2}^{2}|}\big{(}\frac{1}{\beta_{0}}+\ln
t_{0}\big{)}.$

Hence, there exists a number $T_{1}=T_{1}(\xi^{0},\beta_{0})>1$ (see (3.8)),
such, that for all $t_{0}\geq T_{1}$ the function
$sign(\lambda_{1}^{2}-\lambda_{2}^{2})\Re(e^{\prime}_{2,\lambda}(\rho))$ is
negative.

Lemma 2.4 and therefore Theorem 2.5 are proved.

##  4\. An example

To illustrate the theorems proved above consider the following Cauchy problem
(see [30])

(4.1) $\displaystyle D_{\ast}^{\beta_{1}}u_{1}(t,x)$
$\displaystyle=-D^{2}u_{1}(t,x)-Du_{2}(t,x),\quad t>0,\ -\infty<x<\infty,$
(4.2) $\displaystyle D_{\ast}^{\beta_{2}}u_{2}(t,x)$
$\displaystyle=-Du_{1}(t,x)-D^{2}u_{2}(t,x),\quad t>0,\ -\infty<x<\infty,$
(4.3) $\displaystyle u_{1}(0,x)$ $\displaystyle=\varphi_{1}(x),\quad
u_{2}(0,x)=\varphi_{2}(x),\quad-\infty<x<\infty.$

It is not hard to see that the symbol of the operator on the right hand side
of (4.1)-(4.2) is symmetric and has the representation

(4.4) $\mathcal{A}(\xi)=\begin{bmatrix}-\xi^{2}&-\xi\\\
-\xi&-\xi^{2}\end{bmatrix}=\begin{bmatrix}1/2&1/2\\\
-1/2&1/2\end{bmatrix}\begin{bmatrix}-\xi^{2}+\xi&0\\\
0&-\xi^{2}-\xi\end{bmatrix}\begin{bmatrix}1&-1\\\ 1&1\end{bmatrix}.$

As is seen from (4.4) that $\lambda_{1}(\xi)=-\xi^{2}+\xi$ and
$\lambda_{2}(\xi)=-\xi^{2}-\xi.$ The solution $U(t,x)=\langle
u_{1}(t,x),u_{2}(t,x)\rangle$ to Cauchy problem (4.1)-(4.3) has the
representation:

$\displaystyle u_{1}(t,x)$
$\displaystyle=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\left[\frac{1}{2}E_{\beta_{1}}((-\xi^{2}+\xi)t^{\beta_{1}})+\frac{1}{2}E_{\beta_{2}}((-\xi^{2}-\xi)t^{\beta_{2}})\right]\hat{\varphi}_{1}](\xi)d\xi$
$\displaystyle+\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\left[\frac{1}{2}E_{\beta_{1}}((-\xi^{2}+\xi)t^{\beta_{1}})-\frac{1}{2}E_{\beta_{2}}((-\xi^{2}-\xi)t^{\beta_{2}})\right]\hat{\varphi}_{2}](\xi)d\xi;$
$\displaystyle u_{2}(t,x)$
$\displaystyle=\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\left[\frac{1}{2}E_{\beta_{1}}((-\xi^{2}+\xi)t^{\beta_{1}})-\frac{1}{2}E_{\beta_{2}}((-\xi^{2}-\xi)t^{\beta_{2}})\right]\hat{\varphi}_{1}(\xi)d\xi$
$\displaystyle+\frac{1}{2\pi}\int\limits_{-\infty}^{\infty}\left[\frac{1}{2}E_{\beta_{1}}((-\xi^{2}+\xi)t^{\beta_{1}})+\frac{1}{2}E_{\beta_{2}}((-\xi^{2}-\xi)t^{\beta_{2}})\right]\hat{\varphi}_{2}(\xi)d\xi.$

Moreover, obviously, $\lambda_{k}(\xi)\leq 0,\ k=1,2,$ for all $\xi$
satisfying the inequality $|\xi|\geq 1.$ It is not hard to verify, that

$K_{1,1}(\xi,\hat{\Phi}(\xi))=K_{2,1}(\xi,\hat{\Phi}(\xi))=\frac{1}{2}\hat{\varphi}_{1}(\xi)+\frac{1}{2}\hat{\varphi}_{2}(\xi)$

and

$K_{1,2}(\xi,\hat{\Phi}(\xi))=\frac{1}{2}\hat{\varphi}_{1}(\xi)-\frac{1}{2}\hat{\varphi}_{2}(\xi),\quad
K_{2,2}(\xi,\hat{\Phi}(\xi))=-\frac{1}{2}\hat{\varphi}_{1}(\xi)+\frac{1}{2}\hat{\varphi}_{2}(\xi).$

Therefore for the corresponding determinant one has

$\big{|}K_{j,l}\big{(}\xi^{0},\widehat{\Phi}(\xi^{0})\big{)}\big{|}=\frac{1}{2}\big{(}\hat{\varphi}^{2}_{2}(\xi^{0})-\hat{\varphi}^{2}_{1}(\xi^{0})\big{)}$

and condition (2.11) has the form

$\hat{\varphi}^{2}_{2}(\xi^{0})\neq\hat{\varphi}^{2}_{1}(\xi^{0}),\quad\text{or}\quad|\hat{\varphi}_{2}(\xi^{0})|\neq|\hat{\varphi}_{1}(\xi^{0})|.$

In this case the unknown orders $\beta_{1}$ and $\beta_{2}$ are the unique
roots of the following equations

$E_{\beta_{1}}\big{(}\lambda_{1}(\xi^{0})t^{\beta_{1}}_{0}\big{)}=-\frac{d_{1}+d_{2}}{\hat{\varphi}_{1}(\xi^{0})+\hat{\varphi}_{2}(\xi^{0})},$
$E_{\beta_{2}}\big{(}\lambda_{2}(\xi^{0})t^{\beta_{2}}_{0}\big{)}=\frac{d_{1}-d_{2}}{\hat{\varphi}_{2}(\xi^{0})-\hat{\varphi}_{1}(\xi^{0})},$

respectively.

## References

  * [1] Sh. Alimov, R. Ashurov, Inverse problem of determining an order of the Caputo time-fractional derivative for a subdiffusion equation, J. Inverse Ill-Posed Probl. 28 (2020), Issue 5, pp. 651-658.
  * [2] R. Ashurov, Yu. Fayziev, Determination of fractional order and source term in a fractional subdiffusion equation, arXiv:submit/3264960[math. AP]8 Jul 2020.
  * [3] R. Ashurov, S. Umarov, Determination of the order of fractional derivative for subdiffusion equations, Fractional Calculus and Applied Analysis, 2020, 12, pp. 1-17.
  * [4] R.R. Ashurov, R. Zunnunov, Initial-boundary value and inverse problems for subdiffusion equation in $R^{N}$. Fractional Differential Calculus, 2020, 10 (2), 291–306.
  * [5] H. Bateman, Higher transcendental functions, McGraw-Hill, 1953.
  * [6] D. Benson, M. Meerschaert, J. Revielle, Fractional calculus in hydrology modeling: A numerical perspective. Advances in water resources, 51 479–497, 2013.
  * [7] J. Cheng, J. Nakagawa, M. Yamamoto, T. Yamazaki, Uniqueness in an inverse problem for a one-dimensional fractional diffusion equation, Inverse Prob., 4 (2009), 1–25.
  * [8] S. Das, P.K. Gupta, A mathematical model on fractional Lotka-Volterra equations. Journal of theoretical biology, 277 (1), 1-6, 2011.
  * [9] M. M. Dzherbashian, _Integral Transforms and Representation of Functions in the Complex Domain_ (in Russian), M. NAUKA, 1966.
  * [10] Ch. Guo, Sh. Fang, Stability and approximate analytic solutions of the fractional Lotka-Volterra equations for three competitors. Advanced difference equations, 219, 1-14, 2016.
  * [11] Handbook of Fractional Calculus with Applications. Volume 2: Fractional Differential Equations. Editors: Kochubey A., Luchko Yu. De Gruyter, 2019.
  * [12] R. Hilfer, Applications of Fractional Calculus in Physics. World Scientific, 2000.
  * [13] R. Islam, A. Pease, D. Medina, T. Oraby, Integer Versus Fractional Order SEIR Deterministic and Stochastic Models of Measles. International Journal of Environmental Research and Public Health, 17 (6), 1-19, 2020.
  * [14] J. Janno, Determination of the order of fractional derivative and a kernel in an inverse problem for a generalized time-fractional diffusion equation, Electronic J. Differential Equations. 216 (2016), 1-28.
  * [15] N.A. Khan, O. Razzaq, S.P. Mondal, Q. Rubbab Fractional order ecological system for complexities of interacting species with harvesting threshold in imprecise environment. Advances in Difference Equations, 405, 1-34, 2019.
  * [16] A.A. Kilbas, H.M. Srivastava, J.J. Trijillo, Theory and Applications of Fractional Differential Equations. Elsevier Science, 2006.
  * [17] Z. Li, Y. Liu, M. Yamamoto, Inverse problems of determining parameters of the fractional partial differential equations, Handbook of fractional calculus with applications, DeGruyter, 2 (2019), 431–442.
  * [18] Z. Li, Y. Luchko, M. Yamamoto, Analyticity of solutions to a distributed order time-fractional diffusion equation and its application to an inverse problem, Comput. Math. Appl. 73 (2017), 1041-1052.
  * [19] Z. Li, M. Yamamoto, Uniqueness for inverse problems of determining orders of multi-term time-fractional derivatives of diffusion equation, Appl. Anal., 94 (2015), 570–579.
  * [20] T. Machado, A. Lopes, Relative fractional dynamics of stock markets. Nonlinear dynamics, 86 (3), 1613–1619, 2016.
  * [21] R. Magin, Fractional Calculus in bioengineering. Critical reviews in biomedical engineering, 32 (1), 1-104, 2004.
  * [22] F. Mainardi, Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models. Imperial College Press. 2010.
  * [23] R. Metzler, J. Klafter, The random walk’s guide to anomalous diffusion: a fractional dynamics approach. Phys. Rep. 339, no. 1 1–77, 2000.
  * [24] I. Podlubny, Fractional Differential Equations. Academic Press, 1998.
  * [25] Rana S., Bhattacharya S., Pal J., Guerekata G., Chattopadhyay. Paradox of enrichment: A fractional differential approach with memory. Physica A: Statistical Mechanics and its Applications, 392 (17), 3610–3621, 2013.
  * [26] F. Rihan, Numerical Modeling of Fractional-Order Biological Systems. Abstract and Applied Analysis, 2013, 1–13, 2013.
  * [27] S.G. Samko, A.A. Kilbas, O.I. Marichev, Fractional Integrals and Derivatives: Theory and Applications. Gordon and Breach Science Publishers, 1993.
  * [28] E. Scalas, R. Gorenflo, F. Mainardi, Fractional calculus and continuous-time finance. Physica A: Statistical mechanics and its applications, 284 (1-4), 376–384, 2000.
  * [29] S. Umarov, Introduction to Fractional and Pseudo-Differential Equations with Singular Symbols. Springer, 2015.
  * [30] S. Umarov, R. Ashurov, Y. Chen, On a method of solution of systems of fractional pseudo-differential equations. Fractional Calculus and Applied Analysis, to apear.
  * [31] S. Umarov, M. Hahn, K. Kobayashi, Beyond the Triangle: Brownian Motion, Ito Calculus and Fokker-Planck equations - Fractional Generalizations. World Scientific, 2018.
  * [32] B. West, Physics of Fractal Operators. Springer, 2003.
  * [33] A. Zeb, G. Zaman, M.I. Chohan, Sh. Momani, V.S. Erturk, Analytic numeric solution for SIRC epidemic model in fractional order. Asian J. of Math and Appl. 2013, 1-19, 2013.