text:
# The internal structure and formation of early-type galaxies: the gravitational–lens system MG2016+112 at 𝑧=1.004

## 1 Introduction

Understanding the internal structure of galaxies, is crucial to explain the variety of their observed morphologies and the physical processes that created them. It also provides an opportunity to test physics on galactic scales.

Because early-type (E/S0) galaxies are relatively simple systems, dynamical modeling of their internal structure can be attempted with hope for success (e.g. de Zeeuw & Franx 1991; Bertin & Stiavelli 1993; Merritt 1999). Knowledge of this internal structure provides a vital discriminant to the processes involved in their formation (e.g. van Albada 1982), which is currently a highly controversial issue (e.g. Schade et al. 1999).

In particular, in the cold dark matter (CDM) cosmological model, galaxies are embedded in dark matter halos, which dominate the total mass of the system and have a characteristic mass distribution (Navarro, Frenk & White 1997, hereafter NFW; Moore et al. 1998; Bullock et al. 2001; Ghigna et al. 2000). Measuring the distribution of dark matter in E/S0 galaxies would therefore not only shed light on their formation history, but also provide a powerful test of the CDM cosmological model. This test is similar, but complementary in the scales and physical conditions probed, to the still controversial tests performed using the kinematics of low surface-brightness and dwarf galaxies (e.g. McGaugh & de Blok 1998; van den Bosch et al. 2000; Swaters, Madore & Trewhella 2000; Salucci & Burkert 2000; Borriello & Salucci 2001; van den Bosch & Swaters 2001; de Blok et al. 2001; de Blok & Bosma 2002; Jimenez, Verde & Peng 2002).

Unfortunately, E/S0 galaxies generally lack simple dynamical tracers at large radii, such as HI in spiral galaxies, and little is known about their kinematics beyond the region dominated by stellar mass (1–2 effective radii, Re). Hence, the amount and distribution of dark matter within E/S0 galaxies remain to date poorly constrained even in the local Universe (e. g. Bertin et al. 1994; Franx, van Gorkom & de Zeeuw 1994; Carollo et al. 1995; Rix et al. 1997; Gerhard et al. 2001). In a few exceptional cases kinematic tracers at large radii have been found and used to investigate the presence and distribution of dark matter halos (e. g. Mould et al. 1990; Franx et al. 1994; Hui et al. 1995; Arnaboldi et al. 1996). However, these methods do not seem likely to be applicable at significant look-back times and thus to the investigation of the evolution of the internal structure of E/S0 galaxies.

An independent measure of the mass distribution of E/S0 galaxies is provided by gravitational lensing, which is not sensitive to its dynamical state. For example, gravitational lensing alone has been used to show the existence of dark matter in individual E/S0 galaxies (e. g. Kochanek 1995) and recently to constrain the dark matter halo profile from gravitational-lens statistics (e.g. Keeton 2001; Kochanek & White 2001), and galaxy-galaxy lensing signal (e.g. Seljak 2002). The combination of dynamical and lensing analyses can significantly reduce the degeneracies inherent to each method and result in tighter constraints on the dark-matter distribution, as well as on the luminous and total mass distribution. In addition, most lenses are found at intermediate redshift ($`0.1<z<1`$) and therefore offer the possibility to probe the mass distribution at significant look-back times.

### 1.1 The Lenses Structure and Dynamics (LSD) Survey

To investigate the internal structure of E/S0 galaxies and test CDM predictions, we have started a program to measure the kinematic profiles of a sample of E/S0 galaxies that are also gravitational lenses; the Lenses Structure and Dynamics Survey (hereafter the LSD Survey 222see ‘www.astro.caltech.edu/$``$tt/LSD’ or ‘www.its.caltech.edu/$``$koopmans/LSD’). By means of deep spectroscopy, using the Echelle Spectrograph and Imager (ESI) at the W. M. Keck–II Telescope, we are deriving extended kinematic profiles (beyond the effective radius) for the nearest and brightest lenses, along major and sometimes minor axes, while for the most distant and faint ones we are measuring a luminosity-weighted velocity dispersion. The ultimate aim of the LSD Survey is to produce a database of measurements of internal kinematics for the largest possible number of E/S0 lens galaxies, covering a range of redshifts and masses.

The target lenses have to: (i) have E/S0 morphology, as determined from Hubble Space Telescope (HST) images; (ii) be relatively undisturbed by the light of the lensed object at ground based resolution; (iii) be relatively isolated, since a nearby cluster, large group, or massive companion, would introduce additional uncertainties in the measurement of the mass of the lens; (iv) be bright enough that accurate internal kinematics can be determined in a few hours of integration. A preference is given to objects where the redshift of the background source is known. If the redshift of the source is not known and the source is bright enough that the redshift can be determined in a few hours of integration time, the redshift is also determined.

In addition to help clarifying the origin of E/S0 galaxies, other important results will be obtained if a better understanding of the internal structure of E/S0 and its evolution with redshift is achieved. For example we can mention the determination of the Hubble constant from time-delay, and the calculations of lens cross-sections.

### 1.2 A case study: MG2016+112

In a previous paper (Koopmans & Treu 2002; hereafter KT02) we described the Keck and HST observations of the E/S0 lens galaxy D in MG2016+112, the most distant spectroscopically confirmed lens galaxy known ($`z=1.004`$). A review of the current observations of this system can be found in Koopmans et al. (2002; hereafter K02), who also propose a detailed lens model, where the image separation is caused predominantly by galaxy D, and find a central velocity dispersion $`\sigma `$=320–340 km s-1 for an isothermal total mass distribution. This value should correspond closely to the central stellar velocity dispersion (e. g. Kochanek 1994; Kochanek et al. 2000). Our direct measurement ($`\sigma =328\pm 32`$ km s-1; KT02) shows that D is indeed a massive old E/S0 galaxy, possibly embedded in a proto-cluster, with a central stellar velocity dispersion in agreement with that inferred from the lens models.

In this paper we present a detailed analysis of the internal structure of galaxy D, based on the observed velocity dispersion and the exquisite HST images available from the CASTLES database. In particular, we focus on the density profiles of (i) the dark matter halo and (ii) the total luminous plus dark mass distribution. Apart for its high redshift (corresponding to $`8`$ Gyr look-back time333throughout the paper, we assume the Hubble Constant, the matter density, and the cosmological constant to be H0=65 km s-1 Mpc-1, $`\mathrm{\Omega }_\mathrm{m}=0.3`$ and $`\mathrm{\Omega }_\mathrm{\Lambda }=0.7`$, respectively.), this system is particularly interesting because its Einstein radius (REinst=13.7$`\pm `$0.1 kpc) is much larger than the effective radius of the lens galaxy ($`R_\mathrm{e}=2.7\pm 0.5`$ kpc). For this reason, we can use gravitational lensing to probe the mass distribution to much larger radii than is typically possible, even in local E/S0 galaxies.

The paper is organized as follows. First, the luminous and dark-matter models are described in Sect. 2. In Sect. 3, constraints are placed on the dark and luminous matter distributions inside the Einstein radius, and compared to CDM predictions in Sect. 4. In Sect. 5, we use the results found for MG2016+112 to predict observables for other similar systems. Conclusions and a discussion are given in Sect. 6. Throughout this paper, $`r`$ is the radial coordinate in 3-D space, while $`R`$ is the radial coordinate in 2-D projected space.

## 2 Models

Based on the lens models in K02, one finds a total mass enclosed within the Einstein radius, $`M(<R_{\mathrm{Einst}})M_E=1.1\times 10^{12}`$ M. The enclosed mass has been corrected for the slight ellipticity of the lens potential and, because the source lies almost on the observer-lens axis, is nearly independent of the mass distribution assumed for the lens model (e.g. Schneider, Ehlers & Falco 1992; Chapter 8). An error of $`10\%`$ is assumed for the enclosed mass (K02). This error is larger than what one might expect from a range of lens models (e.g. see Kochanek 1991), but is chosen to include any remaining uncertainties related to possibly undetected mass contributions from nearby galaxies. Note, however, that to first order, the contribution from the most massive nearby structures has been corrected for in the lens models (K02). The error estimate on the enclosed mass can thus be regarded as conservative. Later, we will estimate the effect of this potential systematic error on our conclusions and show it to be negligible.

The total luminosity of the system is $`(1.6\pm 0.2)10^{11}`$ L⊙,B (KT02), implying an average mass-to-light ratio M/L$`{}_{B}{}^{}=8\pm 1`$ M/L⊙,B inside the Einstein radius, assuming an $`R^{1/4}`$surface-brightness profile. This value is similar to the values found for the stellar populations of local E/S0 galaxies (e.g. Gerhard et al. 2001). However, considering the much smaller age of the stellar populations of galaxy D – estimated to be $`3.0\pm 0.8`$ Gyr and not larger than the age of the Universe at $`z=1.004`$ ($`6`$ Gyr; KT02) – the stellar M/L ratio of galaxy D must be smaller than in the local Universe, due to luminosity evolution of the stellar population. Given the large value of M/L, we conclude that a significant amount of dark matter is likely to be present inside $`R_{\mathrm{Einst}}`$.

We choose to describe the mass distribution of MG2016+112 with two spherical components, one for the luminous matter and one for the dark matter halo. In particular, we model the luminous mass distribution with a Hernquist (1990) model

$$\rho _L(r)=\frac{M_{}r_{}}{2\pi r(r+r_{})^3},$$
(1)
that well reproduces the $`R^{1/4}`$surface brightness profile for $`r_{}=R_e/1.8153`$ ($`M_{}`$ is the total stellar mass) and is simple to treat analytically (e.g. Ciotti, Lanzoni & Renzini 1996). The dark-matter distribution is modeled as

$$\rho _d(r)=\frac{\rho _{d,0}}{(r/r_b)^\gamma (1+(r/r_b)^2)^{\frac{3\gamma }{2}}}$$
(2)
which closely describes a NFW profile for $`\gamma =1`$, and has the typical asymptotic behavior at large radii found from numerical simulations of dark matter halos $`r^3`$ (e. g. Ghigna et al. 2000). Note that, for $`r_b`$ significantly exceeding $`R_{\mathrm{Einst}}`$, the particular choice of the outer slope has only a negligible effect on the kinematics of the lens galaxy. In particular, for $`r_b\mathrm{}`$ and $`\gamma =2`$, Eq. 2 reduces to a singular isothermal sphere, a functional form often used to describe dark matter halos. In addition, Eq. 2 provides a simple parametrization of the controversial inner slope for $`r<r_b`$, the region better constrained by our observations. The length scale ($`r_b`$) and the density scale ($`\rho _{d,0}`$) determine the virial mass of the dark matter halo (e. g. Bullock et al. 2001).

The velocity dispersion is computed solving the spherical Jeans equation (e.g. Binney & Tremaine 1987), assuming an Osipkov-Merritt (Osipkov 1979; Merritt 1985a,b) parametrization of the anisotropy $`\beta `$:

$$\beta (r)=1\frac{\sigma _\theta ^2}{\sigma _r^2}=\frac{r^2}{r^2+r_i^2},$$
(3)
where $`\sigma _\theta `$ and $`\sigma _r`$ are the tangential and radial component of the velocity dispersion. Note that $`\beta >0`$ by definition and therefore tangentially anisotropic models are not considered in this analysis. Finally, the velocity dispersion as observed through our spectroscopic aperture ($`\sigma _{ap}`$; effectively corresponding to a circular aperture of radius $`0\stackrel{}{\mathrm{.}}652R_e`$) is computed as the average of the projected velocity dispersion weighted by the surface brightness.

## 3 Luminous and dark matter in MG2016+112

### 3.1 Constraints from kinematics, lensing geometry, and surface brightness distribution

The effective radius and the total mass ($`M_E`$) within the Einstein radius are fixed by the observations (KT02). This leaves four free parameters in the model: the inner slope ($`\gamma `$) of the dark matter halo, the length scale of the dark matter component ($`r_b`$), the mass-to-light ratio of the luminous component (M$`{}_{}{}^{}/L_B`$) and the anisotropy radius ($`r_i`$). In Fig. 3.1 we show the likelihood contours444Gaussian error distributions are assumed. of $`\gamma `$ and M$`{}_{}{}^{}/L_B`$, given the observed velocity dispersion ($`\sigma _{ap}=304\pm 27`$ km s-1) for a range of values of $`r_b`$. Since changing $`r_i`$ has very little effect on the likelihood contours, only the ones obtained for $`r_i=R_e`$ are shown (for more discussions see also also the rest of this section and in Sect.5).       Constraints on the inner slope of the dark matter halo ($`\gamma `$) and the stellar mass-to-light ratio ($`M_{}/L_B`$, in solar units). Thin contours represent the 68%, 95%, 99% confidence limits given ME, Re and $`\sigma _{ap}`$. The thick contours represent the same limits, if the independent measurement of $`M_{}/L_B`$ from the evolution of the FP is included (Eq. 6). The subpanels present the results for different values of the dark-matter length scale ($`r_b`$). The contours change negligibly for different values of the anisotropy radius anisotropy radius ($`r_i/R_e=0.5,1,2`$), thus only the case $`r_i/R_e=1`$ is shown. See Sect. 3 for a detailed discussion.  Two main conclusions can be drawn from Fig. 1. First, since M$`{}_{}{}^{}/L_B4`$ M/L⊙,B(99% CL), i.e. smaller than the total $`M/L_B`$=$`8\pm 1`$ M/L⊙,B, dark matter is present within $`R_{\mathrm{Einst}}`$ and contributes $`50\%`$ (99% CL) of the mass. As an independent check, if we assume that $`M_E`$ arises from the luminous component only with a constant M$`{}_{}{}^{}/L_B=8`$M/L⊙,B, we find $`\sigma _{ap}=530`$ km s-1 (virtually independent of the isotropy radius), inconsistent with the observed value at better than 8$`\sigma `$. Second, we find that for large values of $`r_b/R_{\mathrm{Einst}}`$ (e. g. $`r_b/R_{\mathrm{Einst}}=2`$ shown in Fig. 1) and small values of $`M_{}/L_B`$, $`\gamma `$ approaches the value 2. Indeed, results from calculations with M$`{}_{}{}^{}/L_B0`$ and $`r_b\mathrm{}`$, show that $`\gamma =2.0\pm 0.1`$ (1–$`\sigma `$ error). In other words, the total mass can be very well described by a single power-law density distribution $`\rho _tr^\gamma ^{}`$ with effective555We use the term effective slope to emphasize that this is the slope that we get with a single power law. The real mass distribution could change slope between the constraints given within the spectroscopic aperture radius and at $`R_{Einst}`$, as long as the effective trend is preserved. slope $`\gamma ^{}=2.0\pm 0.1`$. A systematic error of 10% on $`M_\mathrm{E}`$ changes the inferred slope by only 0.05, whereas a 20% error (1–$`\sigma `$) on $`R_\mathrm{e}`$ affects the result by only 0.02. Similarly, a 20% error (1–$`\sigma `$) on the assumed aperture radius of $`0\stackrel{}{\mathrm{.}}65`$ changes the result by 0.03. Even a factor 10 deviation from $`r_i=R_\mathrm{e}`$ results in $`0.05`$ error, although the effect on the central velocity dispersion666defined as the luminosity-weighted velocity dispersion within a circular aperture of radius $`R_\mathrm{e}/8`$ $`\sigma `$ in that case is much more profound, with significantly larger values of $`\sigma `$ for smaller $`r_i`$ (i.e. more radial orbits, see Fig. 4). Hence, $`\gamma ^{}`$ is robust against reasonable changes in the parameters $`r_i`$, $`R_\mathrm{e}`$, $`M_\mathrm{E}`$, the aperture radius, and $`r_bR_{\mathrm{Einst}}`$.

We can make a very simple argument to understand why the effective slope is so tightly constrained by our observations. For a total mass profile $`\rho _\gamma ^{}r^\gamma ^{}`$ the enclosed mass is $`M(<r)r^{3\gamma ^{}}`$. Therefore, the mass within the effective radius ($`M_e`$) depends on the mass within the Einstein radius ($`M_E`$) as

$$M_e=M_E\left(\frac{R_e}{R_{Einst}}\right)^{3\gamma ^{}}M_Ex_E^{\gamma ^{}3},$$
(4)
where $`x_ER_{Einst}/R_e`$ has been introduced for simplicity of notation. Hence, assuming $`M_e\sigma ^2`$ for a given $`R_e`$,

$$\delta \gamma ^{}\frac{2\delta \sigma }{\sigma |\mathrm{log}(x_E)|},$$
(5)
where $`\delta \gamma ^{}`$ and $`\delta \sigma `$ are the errors on the effective slope and luminosity-weighted velocity dispersion, respectively. Using the values for MG2016+112 in Eq. 5 we recover $`\delta \gamma ^{}0.1`$ in good agreement with the error inferred from our model.

### 3.2 Further constraints: the evolution of the Fundamental Plane

Recent studies have shown that E/S0 galaxies both in clusters and in the field define a tight Fundamental Plane (FP; Dressler et al. 1987; Djorgovski & Davis 1987) out to $`z0.70.8`$ (van Dokkum et al. 1998; Treu et al. 1999, 2002; see also Kochanek et al. 2000), with slopes very similar to the ones observed in the local Universe. Assuming that galaxy D in MG2016+112 lies on a FP with slopes as in the local Universe, we can obtain an additional constraint on its internal structure. In fact, the evolution of the intercept of the FP with redshift can be related to the evolution of the average effective mass to light ratio ($`\mathrm{\Delta }\mathrm{log}M/L_B`$; see e. g. Treu et al. 2001 for discussion). For D we obtained $`\mathrm{\Delta }\mathrm{log}M/L_B=0.62\pm 0.08`$, i.e. intermediate between the cluster and field value (KT02). We can use this measurement to infer $`M_{}/L_B`$ of D assuming that

$$\mathrm{log}(M_{}/L_B)_z=\mathrm{log}(M_{}/L_B)_0+\mathrm{\Delta }\mathrm{log}(M/L_B),$$
(6)
where the second term on the right hand side of the equation is measured from the evolution of the FP. The first term on the right hand side of the equation can be measured for local E/S0 galaxies: for example using the data from Gerhard et al. (2001; and references therein), we find a good correlation between $`\sigma ^2R_e`$ and $`M_{}/L_B`$. Using the values of MG2016+112, this correlation yields $`M_{}/L_B(z=0)=7.3\pm 2.1`$ M/L⊙,B, where the error includes the contribution from the scatter between dynamical and stellar population measurements777A very similar result is obtained by taking simply the average stellar M/L of the sample $`7.8\pm 2.7`$M/L⊙,B.. Hence, from Eq. 6 we find $`M_{}/L_B(z=1.004)=1.8\pm 0.7`$ M/L⊙,B. We note two things: (i) the value $`M_{}/L_B=8`$ M/L⊙,B required to explain the enclosed mass solely by luminous matter is inconsistent with the above, independently–derived, value at the 9–$`\sigma `$ level, and (ii) the allowed range of $`M_{}/L_B`$ found from our dynamical modeling (Fig. 1) overlaps very nicely with the value derived from the evolution of the FP. Using this additional constraint we obtain the probability contours plotted as thick solid lines in Fig. 3.1. Again, the results depend very modestly on the anisotropy radius, while the slope $`\gamma `$ becomes steeper as $`r_b`$ increases. In general, $`\gamma =2`$ can be ruled out at better than 2–4 $`\sigma `$ (depending on $`r_b`$), and therefore the dark matter profile is not isothermal.       Dark and luminous matter best-fit density (in units of $`10^{10}`$ M$`_{}`$kpc-3) and velocity dispersion profiles (in km s-1) as function of radius (in kpc). Lower panel: Density profile of the luminous component (dot-dashed), of the dark matter halo (dotted) and the sum of the two (solid) for several characteristic best-fitting models. An $`r^2`$ density profile (dashed) reproduces the total density profile remarkably accurately. Middle panel: Relative residuals of the total mass distribution ($`\rho _{tot}`$) with respect to power laws $`\rho _\gamma ^{}=r^\gamma ^{}`$, with $`\gamma ^{}=2.2,2,1.8`$ from top to bottom. Upper panel: Luminosity-weighted velocity dispersion profile as a function of circular aperture. The solid lines are obtained for the best-fitting two-component model. The thin dot-dashed lines correspond to power-law mass distributions with $`\gamma ^{}=`$2.2,2,1.8 from top to bottom. The error bar represents our measurement $`\sigma _{ap}=304\pm 27`$km s-1. We assume $`M_{}/L_B=1.8`$M/L⊙,B.

Let us now examine the total mass distribution of the two-component model in more detail with the help of the quantities plotted in Fig. 2. First (lower panel), we show the density profiles of a characteristic best-fitting model ($`r_b=2R_e`$, $`r_i=R_e`$; $`M_{}/L=1.8`$ M/L⊙,B). Even though luminous and dark matter themselves are not isothermally distributed, the total mass distribution is very close to isothermal, confirming what we found in Section 3.1. Second (middle panel), we show the relative difference between a power law mass profile $`\rho _\gamma ^{}r^\gamma ^{}`$ and the total mass profile obtained with our two-component model. A profile with $`\gamma ^{}=2`$ describes the total mass profile very well. By contrast, $`\gamma ^{}=1.8`$ or $`\gamma ^{}=2.2`$ result in significantly worst representations of the total mass profile. Third (upper panel), we compare the luminosity-weighted velocity dispersion profile of our two-component model with those obtained for $`\gamma ^{}=`$1.8, 2.0 and 2.2. The $`\gamma ^{}=2`$ power law mass profile and the two-component model produce almost identical velocity dispersion profiles, while $`\gamma ^{}=1.8,2.2`$ are cleary inconsistent with the observed $`\sigma _{\mathrm{ap}}`$ (point with error bar in the plot). For completeness, we note that we find $`\sigma /\sigma _{ap}=1.06`$, consistent with our assumption of $`1.08\pm 0.05`$ in KT02.

Finally, we check that our best-fitting models correspond indeed to physical solutions in the sense that they can be generated by a positive distribution function (i.e. they are consistent models, see e.g. Ciotti 1999 and references therein for discussion). For our choice of values of $`r_i/R_e`$, a luminous Hernquist model embedded in a total isothermal mass distribution always satisfies the sufficient condition given in Eq. 6 of Ciotti (1999), and therefore our best-fitting models are consistent.

## 4 Comparison with CDM models

In order to make a proper comparison with CDM predictions, we have to consider that the universal profiles of dark matter halos can be significantly modified by the formation of a galaxy within the halo. A simple treatment of this phenomenon can be given by assuming the so-called “adiabatic approximation” (Blumenthal et al. 1986; see also Mo, Mao & White 1998). Under this approximation, for any given initial mass distribution (baryons are assumed to follow dark matter at the beginning) we can calculate the final dark matter distribution after the galaxy has collapsed, for example to an Hernquist mass distribution (e. g. Keeton 2001; see also Rix et al. 1997). Here, we use the approximation backwards to recover the initial dark matter density profile. In order to facilitate comparison to CDM predictions, we fit the same functional form as in Eq. 2 to the initial density profile by minimizing the fractional difference in density integrated over 0–100 kpc, and derive the initial length scale $`r_{b,i}`$ and an initial slope $`\gamma _i`$. Although it was not guaranteed a priori, we find that this form gives a reasonable description of the initial density profile as well.

Before proceeding to the calculation, we reduce the degrees of freedom of the problem by noticing that CDM models indicate that the length scale of systems as massive as galaxy D should be significantly larger than $`15`$ kpc, even at $`z1`$ (Bullock et al. 2001). Therefore, in the spirit of the comparison with CDM predictions, we can limit our parameter space to $`r_b>R_{\mathrm{Einst}}`$, where the results depend very little on the exact value of $`r_b`$. The length scale does not change significantly ($`<30`$% depending on the slope for $`r_{bi}<40`$ kpc) during collapse, therefore our assumption on $`r_b`$ being large because $`r_{b,i}`$ is large is justified in this context. For definiteness, we pick $`r_b=2R_{\mathrm{Einst}}`$ so that the observed profiles are shown in Fig. 3.1 and 3.2. With this assumption, the initial slope $`\gamma _i`$ depends only on the observed slope $`\gamma `$ and on $`M_{}/L_B`$.

In Fig 4 we show the contour levels of $`\gamma _i`$ in the $`\gamma M_{}/L_B`$ plane. As $`M_{}/L_B0`$ , the effect of baryons becomes negligible and $`\gamma \gamma _i`$. On the other hand, as $`M_{}/L_B`$ increases, the collapse becomes increasingly important and flatter initial dark matter halos correspond to increasingly steeper final halos. To constrain $`\gamma _i`$ with the observations, we overlay the likelihood contours of $`\gamma M_{}/L_B`$ derived from our dynamical model (Fig. 1), and consider only the portion of the $`\gamma _i`$ contours enclosed by the likelihood contours. Considering only monotonically decreasing profiles, the range $`\gamma _i=0.01.4`$ maps onto a very tight range of final slopes $`\gamma 1.301.75`$. In conclusion – within the context of this simple adiabatic contraction model – $`\gamma _i=1.5`$ is only marginally consistent with our observations, whereas $`\gamma _i=1`$ fits within the 68% likelihood contour.       Relation between the final ($`\gamma `$) and the initial ($`\gamma _i`$) inner slope of the CDM dark-matter halo of galaxy D in MG2016+112, and the stellar mass-to-light ratio ($`M_{}/L_B`$), based on a simple adiabatic contraction model (see text). The loci corresponding to different initial slopes are shown, labeled by $`\gamma _i`$. In particular, the two thick contours labeled 1.0 and 1.5 (=$`\gamma _i`$) indicate the NFW and Moore et al. (1998) profiles, respectively. The probability contours given the observations of MG2016+112 for a model with $`r_\mathrm{d}=2R_{\mathrm{Einst}}`$ and $`r_i=R_\mathrm{e}`$ (see Fig.1) are overplotted and labeled by the probability level. As in Fig. 1 solid lines include the constraint from the FP evolution, dashed lines do not.

## 5 Kinematics of Hernquist luminosity density profiles in a logarithmic gravitational potential

In Sect. 3 we have shown that the total mass distribution of galaxy D is well described by an isothermal sphere in the range $``$ 1–15 kpc. Here, we generalize the results found for galaxy D to compute the velocity dispersion profiles that we expect to find for other lens E/S0 galaxies.

Let us assume that E/S0 galaxies have total mass profiles that follow $`\rho _{tot}r^2`$ between a fraction of an effective radius and a few effective radii. Solving the Jeans Equation for a Hernquist luminosity-density profile embedded in a logarithmic potential, we can compute the velocity dispersion profile, depending only on the effective radius ($`R_e`$), the “Einstein velocity dispersion” ($`\sigma _{\mathrm{Einst}}\sqrt{GM_\mathrm{E}/R_{\mathrm{Einst}}})`$, and the anisotropy radius ($`r_i`$) .

In Fig. 5, we show the luminosity-weighted velocity dispersion profile (in units of $`R_e`$ and $`\sigma _{\mathrm{Einst}}`$), for several values of $`r_i`$. Note that the profiles approach the same curve for large $`R/R_e`$ independently of $`r_i`$. Hence, the effect of changing $`r_i`$ on the luminosity-weighted velocity dispersion beyond the effective radius, is not as dramatic as in the inner portions of E/S0 galaxies. For this very reason, the luminosity-weighted velocity dispersion does not firmly constrain $`r_i`$. On the other hand, the velocity dispersion profile itself is much more sensitive to $`r_i`$. Therefore, a spatially resolved measurement of $`\sigma (r)`$, extending from a fraction of to a few $`R_e`$, will allow us to constrain both $`r_i`$ and the total mass distribution.       Luminosity-weighted velocity dispersion profile for a Hernquist luminosity-density profile in logarithmic potential. The axes are scaled in units of $`R_e`$ and $`\sigma _{\mathrm{Einst}}`$ (see text). From top to bottom curves are shown for $`r_i/R_e`$=$`0.2`$,$`0.5`$,1,2,$`\mathrm{}`$ (where $`\mathrm{}`$ corresponds to isotropic models). Note the smaller effect of changing anisotropy radius $`r_i`$, as compared to changing the effective slope $`\gamma ^{}`$ of the total mass distribution (upper panels in Fig. 3.2; to facilitate comparison the two figures show the same range of velocity dispersion and radius). Our measurement of MG2016+112 is shown as a filled circle. The value for 0047-0281 using data from Warren et al. (1998) and Kochanek et al. (2000) is shown as an open triangle.

To our knowledge, the only system (besides MG2016+112) satisfying the LSD requirements888Plus showing no evidence for rotation, since this dynamical model is inappropriate for rotating systems. (Sec 1.1), for which velocity dispersion of the lens galaxy has been measured is 0047-281 (Warren et al. 1996, 1998; the effective radius and image separation are measured by Kochanek et al. 2000). As shown in Figure 5, the observed velocity dispersion of the lens galaxy in 0047-281 agrees very well with the value computed with our model. In future papers we will present velocity dispersion (profiles) for other systems, and further test this model.

## 6 Summary & Discussion

### 6.1 Summary

We have shown that the combination of constraints from stellar kinematics and gravitational lensing is extremely powerful in determining the mass profiles of luminous and dark matter in E/S0 galaxies even at large cosmological distances. In particular, for galaxy D in MG2016+112 ($`z=1.004`$) we have demonstrated that:

(i) Dark matter is unambiguously present and accounts for $`>`$50% (99% CL) of the mass within the Einstein radius. A mass distribution following light with constant $`M/L`$ can be ruled out at more than 8–$`\sigma `$ level. Considering an additional constraint on $`M_{}/L_B`$ from the evolution of the FP, the smallest allowable fraction of dark matter increases to 60% (99% CL).

(ii) The total mass distribution (dark+luminous) inside the Einstein radius is well described by an isothermal sphere, i.e. $`\rho _{tot}r^\gamma ^{}`$ with $`\gamma ^{}=2.0\pm 0.1`$. This result is robust and does not depend significantly on the dynamical state of the system or on the modeling of the mass distribution (one or two component models). The result is only weakly dependent on errors on $`M_E`$, $`r_i`$, $`R_\mathrm{e}`$ and the aperture radius, which in quadrature add at most another 0.1 to the error.

(iii) Using two-component (dark and luminous) spherical models, we find that the inner slope of the dark matter halo $`\gamma <2.0`$ at the 95% CL, i. e. the halo is not isothermal inside the Einstein radius.

(iv) Using a simple adiabatic contraction model, we relate the observed inner slope $`\gamma `$ to the inner slope $`\gamma _i`$ of the dark matter halo before collapse. In particular we find the upper limit $`\gamma _i<1.4(2.0)`$ at 68(99)% CL, marginally consistent with the high resolution CDM simulations that indicate $`\gamma _i`$$``$1.5 (Ghigna et al. 2000; but see Power et al. 2002 for a discussion of the related uncertainties).

These results are robust because the stellar velocity dispersion and gravitational lensing provide two independent mass measurements at two well-separated radii, the effective radius (2.7 kpc) and the Einstein radius (13.7 kpc $`5R_\mathrm{e}`$). This leads to a very well-defined slope inside the Einstein radius, and to stringent limits on the fraction of dark matter.

Higher signal to noise – and possibly spatially resolved observations – are needed to improve these constraints. However, this is unlikely to be feasible with the current generation of space or ground-based telescopes and will probably need the Next Generation Space Telescope or large ground-based telescopes with adaptive optics. Nevertheless, when the measurement and analysis of the dozen systems targeted by the LSD Survey will be completed, we expect to be able to obtain more stringent constraints by combining the results from the individual systems.

Finally, in Sect. 5 we have provided a simple test of whether an isothermal total mass distribution is a good description of other lens galaxies. In particular, we compute the luminosity-weighted velocity dispersion profile for a luminous Hernquist component in a logarithmic potential, which depends very mildly on the anisotropy radius, especially at radii larger than $`R_e`$. The measured velocity dispersion of the only system available besides MG2016+112 (0047+281) agrees very well with the curve computed with this model. As new velocity dispersion measurements of E/S0 lens galaxies become available, they can be readily compared to the expectations of this simple model to test whether indeed the isothermal mass profile is a general feature of the E/S0 population.

### 6.2 Discussion

In this paper and KT02 we presented the following simple picture of the lens galaxy D in MG2016+112. The stellar populations are old and metal rich. A significant amount of dark matter is present and the mass distribution of the luminous and dark component are well constrained, reproducing all the available observations.

However, the good agreement of the total mass profile with an $`r^2`$ power law calls for a physical explanation. This seems not to be a mere coincidence, given that a similar behavior has been observed in some local ellipticals with kinematic tracers at extraordinary large radii (e. g. Mould et al. 1990; Franx et al. 1994; Rix et al. 1997). Whether this galaxy formed by monolithic collapse, or by mergers of smaller subunits in a hierarchical picture, why do these processes conspire to produce an apparent isothermal mass distribution within REinst? In the case of adiabatic contraction, there appears to be no reason why gas could not continue to contract or stop earlier and give rise to a total mass distribution that is either steeper or shallower than isothermal. Thus, formation by adiabatic contraction seems not to be a satisfactory explanation, and it is worth seeking a more fundamental explanation. By contrast, complete violent relaxation produces isothermal distribution functions (Lynden-Bell 1967; Shu 1978), and would explain quite naturally the observations. However, this only applies exactly for an infinite mass distribution. When finite mass distributions are considered, (incomplete) violent relaxation leads to an $`R^{1/4}`$ surface-brightness profile (e.g. van Albada 1982; Hjorth & Madsen 1991, 1995), which goes as $`r^2`$ within the half mass radius, the region of interest here. For the two-component case (luminous and dark), the process should be such as the two components interact to produce the total $`r^2`$ density profile, while preserving the segregation of the luminous component in a centrally concentrated density profile (see e.g. the two-component models of Bertin, Saglia & Stiavelli 1992; a similar outcome can be produced by merging of galaxies with massive halos, see e. g. Barnes 1992; Barnes & Hernquist 1996). In the CDM scenario, it is not clear whether the chaotic changes in the gravitational potential during the collapse of halos satisfy the conditions for violent relaxation (e.g. Flores & Primack 1994 and references therein), although it might not be inconceivable that violent relaxation of the inner parts of E/S0 galaxies is close to complete. A complete discussion of these issues goes beyond the aims of this paper. However, we would like to emphasize that whatever physical processes are involved, the mass distribution of galaxy D in MG2016+112 seems already relaxed $`8`$ Gyr ago.

Finally, we note that our poor knowledge of the slope of the mass profile of lens galaxies is the major systematic uncertainty in the determination of H0 from time-delays (e.g. Koopmans & Fassnacht 1999; Romanowsky & Kochanek 1999; Koopmans 2001). The assumption of isothermal mass distribution leads to $`H_06070`$ km s-1 Mpc-1 (e.g. Koopmans & Fassnacht 1999; Koopmans 2001), in good agreement with other independent measurements. The result for MG2016+112 seems to indicate that the systematic error – the deviation from isothermality – might indeed be very small for E/S0 lens galaxies. However, the LSD Survey, by measuring the average effective slope of the mass profile and its scatter, will allow us to pin down the largest systematic uncertainty in the global determination of H0 from gravitational-lens time delays.

We are grateful to E. Agol, A. Benson, G. Bertin, R. Blandford, L. Ciotti, R. Ellis, M. Stiavelli for the stimulating conversations and insightful comments that helped in shaping this manuscript. Uros Seljak is thanked for useful conversations. We thank the anonymous referee for his prompt and careful report. We acknowledge financial support from NSF and HST grants (AST–9900866; STScI–GO 06543.03–95A; STScI-AR-09222). Finally, the authors wish to recognize and acknowledge the very significant cultural role and reverence that the summit of Mauna Kea has always had within the indigenous Hawaiian community. We are most fortunate to have the opportunity to conduct observations from this mountain.


format: text

text:
# 1 Introduction

## 1 Introduction

In recent years it becomes obvious a reanimation of intense interest to further development of exact methods for solving differential equations. The main reasons here are not only in expansion of science areas where differential equations play prominent role, but also in impetuous development of computer algebra systems (see ).

There are several universal approaches to solving DE. The most promising of them are Lie symmetry methods (see  ) and advanced methods for the linear equations which are based on differential Galois theory  with the first usable algorithm developed by Kovacic , . But in computer algebra very often underlying mathematics is far from the methods normally taught to mathematics students for the same problem .

The aim of this contribution is to show on examples of first-order equations the possibilities for solving differential equations with new method based on Dyson’s operator solution representation in form with time-ordered exponent. This method allows to find instantly the exact operator solution if we succeed in reformulating given problem to the first-order linear one. This solution is formal enough. The next step of the approach is directed to transformation the operator solution to combinations of some operators with relatively simple properties, which would permit to carry out all the operators and obtain as a result the solution of equation ”by quadratures” or to more practical, calculable or usual (standard) expression.

The main advantages of the method are in its compactness, clarity and simpleness. Since we often can obtain some form of exact solution and due to the fact that it is much easier to study the ”prepared frog”. It is essential as well, that the operator method is algorithmic.

It is clear that there are some correlations between operator method and others common methods. But we believe that in introductory paper it is untimely to analyse such correlations.

We choose the first-order equation as a testing area for the operator method because it is most important differential equation in many branches of science and because main aspects of the method are quite easily generalized for other cases.

There are some known solvability conditions for the first-order equations (see for example ). But in the situation, when we are far from be able to solve any equation, the pursuit of new such conditions is very important.

To avoid induced complications we choose a nonrigourous manner of exposition inasmuch we always can verify obtained conclusions by substitution to initial equation. Analytical restrictions are obvious from context as a rule.

The paper is organized as follows. In Sec.2 we introduce basic ideas of the operator method. The applicability of the operator method for non-linear problems is discussed in Sec.3 on example of first-order non-linear ODE. The aim of this section is to introduce on particular case what we will do for the general cases. In Sec.4 we consider an idea of CAS solving algorithm, which is based on operator solution. We believe that this type of direct solving algorithms can play a noticeable part in future CAS developments. The main tools of operator method: the transformation identities are briefly described in Sec.5 and important properties of generalized shift operators are considered in Sec.6. In the Sec.7 we demonstrate algebraical way to simplify operator solutions. Such type of solutions can found more efficient direct CAS algorithms. The expression of the general solution of the equation through particular solution is considered in Sec.8. This section carries a methodical load too. Here we demonstrate the possibilities of operator method in obtaining solutions for more complicated equations. It is shown in Sec.9 that the solution of DE can be expressed in terms of arbitrary function. Among another consequences it brings some recipes for finding of solvability conditions that allow to enlarge DE solving abilities of existent CAS.

## 2 Base of the operator method

The main point of the operator method is based on solution of the *first-order linear* equation obtained by Dyson in the form with time-ordered exponent.

For the *first-order* equation with respect to $`t`$ (in ordinary or partial derivatives) for $`u(t,\stackrel{}{\rho })`$ :

$$\frac{u(t,\stackrel{}{\rho })}{t}=\widehat{L}(t,\stackrel{}{\rho })u(t,\stackrel{}{\rho })+f(t,\stackrel{}{\rho })$$
(1)
with the initial condition

$$u(a,\stackrel{}{\rho })=v(\stackrel{}{\rho }),$$
(2)
where $`\widehat{L}`$ is a *linear operator* (independent variable $`\stackrel{}{\rho }`$ is one dimensional or many dimensional), the solution is given by $`(ta)`$

$$u(t,\stackrel{}{\rho })=\widehat{T}\mathrm{exp}\{_a^t𝑑\tau \widehat{L}(\tau ,\stackrel{}{\rho })\}v(\stackrel{}{\rho })+_a^t𝑑\xi \widehat{T}\mathrm{exp}\{_\xi ^t𝑑\zeta \widehat{L}(\zeta ,\stackrel{}{\rho })\}f(\xi ,\stackrel{}{\rho }).$$
(3)
Here the operator $`\widehat{T}`$ is a ”chronological” (ordering) Dyson’s operator, representing a product, ordered (in order of unincrease) along the co-ordinate $`t`$ , of the operators, which appear in the integrands in the expansion (see below) of the exponent in an operator series:

$`\widehat{T}[\widehat{L}(\tau _1,\stackrel{}{\rho })\mathrm{}\widehat{L}(\tau _m,\stackrel{}{\rho })]=\widehat{L}(\tau _{\alpha _1},`$ $`\stackrel{}{\rho })\widehat{L}(\tau _{\alpha _2},\stackrel{}{\rho })\mathrm{}\widehat{L}(\tau _{\alpha _m},\stackrel{}{\rho }).`$
$`\tau _{\alpha _1}\tau _{\alpha _2}\mathrm{}\tau _{\alpha _m}`$
This chronological operator has to be introduced because the operators in the integrands with different values of the variable $`\tau `$ may not commute. The solution of the problem (1), (2) for $`t<a`$ is similar to (3), but instead of the operator $`\widehat{T}`$ it is necessary to place the operator $`\widehat{T}_0`$ , which represents a product, ordered in reverse order.

We can define chronological operator exponent by another way, e.g., through its differential properties or by way of expansion in an operator series. We can find that

$$\widehat{T}\mathrm{exp}\{_a^t𝑑\tau \widehat{L}(\tau ,\stackrel{}{\rho })\}=\widehat{1}+_a^t𝑑\tau \widehat{L}(\tau ,\stackrel{}{\rho })\widehat{T}\mathrm{exp}\{_a^\tau 𝑑\xi \widehat{L}(\xi ,\stackrel{}{\rho })\},$$

were $`\widehat{1}`$ is the identity operator. Father iterations lead to

$`\widehat{T}\mathrm{exp}\{`$ $`{\displaystyle _a^t}d\tau \widehat{L}(\tau ,\stackrel{}{\rho })\}=\widehat{1}+{\displaystyle }_a^td\tau \widehat{L}(\tau ,\stackrel{}{\rho })+\mathrm{}`$
$`\mathrm{}+{\displaystyle _a^t}𝑑\tau _1{\displaystyle _a^{\tau _1}}𝑑\tau _2\mathrm{}{\displaystyle _a^{\tau _{n1}}}𝑑\tau _n\widehat{L}(\tau _1,\stackrel{}{\rho })\widehat{L}(\tau _2,\stackrel{}{\rho })\mathrm{}\widehat{L}(\tau _n,\stackrel{}{\rho })+\mathrm{}`$ (4)
The most important here is the fact that (3) is really the exact solution of the equation (1) with the initial condition (2). One can very easily verify it by means of differentiation of (3) with respect to $`t`$ .

It is obvious that

$$\widehat{T}\mathrm{exp}\{_a^t𝑑\tau \widehat{L}(\tau ,\stackrel{}{\rho })\}\widehat{T}_0\mathrm{exp}\{_t^a𝑑\tau \widehat{L}(\tau ,\stackrel{}{\rho })\}=\widehat{1}$$

and, e.g., when $`atb`$

$$\widehat{T}\mathrm{exp}\{_a^b𝑑\tau \widehat{L}(\tau ,\stackrel{}{\rho })\}\widehat{T}_0\mathrm{exp}\{_b^t𝑑\tau \widehat{L}(\tau ,\stackrel{}{\rho })\}=\widehat{T}\mathrm{exp}\{_a^t𝑑\tau \widehat{L}(\tau ,\stackrel{}{\rho })\}.$$

It must be stressed here that we do not involve any premises about nature of $`u(t,\stackrel{}{\rho })`$ , it can be both scalar and vectorial function, and it can depend on any number of independent arguments. For operator $`\widehat{L}`$ (and its ”filler”) as well - the only fundamental requirement - it must be a linear operator. Very often non-linear problem is reduced to linear differential equation. If we succeed in reformulating given problem to first-order linear one, we then instantly can find its solution in operator form.

At first sight the chronological operator exponent is a very specific case of linear operators. The following proposition asserts the opposite. Let we have arbitrary but *linear differentiable reversible* operator $`\widehat{A}(t)`$ . From obvious identity

$$\frac{d\widehat{A}(t)}{dt}=[\frac{d\widehat{A}(t)}{dt}\widehat{A}^1(t)]\widehat{A}(t),\widehat{A}(t)|_{t=a}=\widehat{A}(a)$$

we obtain

$$\widehat{A}(t)=\widehat{T}\mathrm{exp}\{_a^t𝑑\tau [\frac{d\widehat{A}(\tau )}{d\tau }\widehat{A}^1(\tau )]\}\widehat{A}(a),$$

so the chronological operator exponent can represent *any* linear differentiable reversible operator, only this form is exactly adjusted for solving defined DE. For any operator of such type there is corresponding linear DE. We can stress here that to solve a DE in explicit form we have to express explicitly the action of corresponding linear operator. It is the challenge: far from all linear operators can be specified by explicit action.

In present paper for simplicity we restrict consideration by first-order DE’s.

## 3 Operator solutions for first-order non-linear ODE

First-order differential equation of arbitrary form

$$\frac{du(t)}{dt}=f(t,u(t))$$
(5)
on account of non-linearity is not suited immediately for solving by operator method. Nevertheless there are many possibilities to convert the problem to linear one. Most of them are concerned with introduction of spaces of larger dimensions. Here we demonstrate some of them, which are almost generally applicable.

If we introduce a new function

$$S(t,\omega )=e^{\omega u(t)},$$

from (5) one can easily derive the following equivalent to (5) first-order (with respect to $`t`$) linear DE for $`S(t,\omega )`$

$$\frac{S(t,\omega )}{t}=\omega f(t,\frac{}{\omega })S(t,\omega ),S(a,\omega )=e^{\omega u(a)},$$
(6)
which formal solution follows immediately from (3) in the form

$$S(t,\omega )=\widehat{T}\mathrm{exp}\{_a^t𝑑\tau \omega f(\tau ,\frac{}{\omega })\}e^{\omega u(a)}.$$
(7)
The solution of equation (5) is obvious from (7):

$$u(t)=\frac{S(t,\omega )}{\omega }|_{\omega =0}.$$

We have to note here that when we look for general or any particular solutions we can utilize indefinite integrals in operator solutions.

The operator solutions of type (3), (7) are formal enough. Meantime, solutions satisfy the following requirements to be the exact solution:

a). Unknown is on the left side. The right side consists of only and only knowns (variables, coefficients, functions, operators and etc.).

b). Direct substitution of the solution reduces the equation to an identity.

The theorem for uniqueness of solution of differential equation shows that if we have an exact solution of a problem in several forms, all the forms can be transformed to each other. So the main problem that remains in this approach is how to transform of the operator solution to that form which can be considered as usable.

The most part of this paper is devoted to hunting such ways.

## 4 An idea of solving algorithm

It would be rather well if we can find a universal algorithm, which would solve any DE. The most attractive feature of the operator method is unusual easiness in obtaining of formal (operator) solutions of DE’s. Can we find on this way a general algorithm to obtain a standard (non-operator) solutions of DEs?

Let us consider *an idea* of such algorithm on special example of solution of equation for unknown function $`u(t)`$:

$$\frac{du(t)}{dt}=u^2(t)+\frac{da(t)}{dt}a^2(t),$$
(8)
where $`a(t)`$ is some known function. Such type of equation is very suitable for *testing* algorithms. First of all since its particular solution is obviously known: $`u(t)=a(t)`$. In the second place since existent CA systems *do not solve* at present most cases of such type of equations (with specific expressions for $`a(t)`$).

The operator solution of equation (8) in form (7) is

$$u(t)=[\frac{}{\omega }\widehat{T}\mathrm{exp}\{𝑑t\omega [\frac{^2}{\omega ^2}+\frac{da(t)}{dt}a^2(t)]\}e^{\omega c}]_{\omega =0},$$
(9)
were $`c`$ is an arbitrary constant (farther on in this section for simplicity let $`c=0`$, i.e. we will seek here only a particular solution of equation). If we represent the chronological exponent as the series (4):

$$u(t)=[\frac{}{\omega }\underset{i=1}{\overset{\mathrm{}}{}}S_i(\omega ,t)]_{\omega =0},$$

terms of series obviously obey the following recurrence relations

$$S_1(\omega ,t)=1,S_i(\omega ,t)=𝑑t\omega [\frac{^2}{\omega ^2}+\frac{da(t)}{dt}a^2(t)]S_{i1}(\omega ,t).$$

Now let us use the following trick. If we compute the sequence of

$$u_n(t)=[\frac{}{\omega }\underset{i=1}{\overset{n}{}}S_i(\omega ,t)]_{\omega =0},$$

we find

$`u_1(t)`$ $`=a(t){\displaystyle 𝑑ta^2(t)},`$
$`u_3(t)`$ $`=a(t)2{\displaystyle 𝑑ta(t)𝑑ta^2(t)}+{\displaystyle 𝑑t[𝑑ta^2(t)]^2},`$
$`u_5(t)`$ $`=a(t)3{\displaystyle 𝑑t[𝑑ta^2(t)]^2}+6{\displaystyle 𝑑t𝑑ta(t)[𝑑ta^2(t)]^2}+\mathrm{}`$
If compare now options in sequence of $`u_n(t)`$ we can find that from some $`n_0`$ an unaltered (survived) part in $`u_n(t)`$ become apparent (the $`a(t)`$ here). It is reasonable that we can suppose that it is just the sought solution. On solution of DEs we always have the opportunity to verify it. If we get affirmative answer, so it is the end of procedure, else we proceed with larger $`n`$.

We have prepared a computer algebra prototype of some variants of the algorithm using the Maple system. The algorithm manages, e.g. the following DE’s (which we choose mostly by chance):

$`DE1:={\displaystyle \frac{du(t)}{dt}}=`$ $`u^2(t)+{\displaystyle \frac{1}{t^{138}}}+41552420t^{90}+4410t^{89}+{\displaystyle \frac{7225}{t^{136}}}13430t^{23}+`$
$`+8330t^{22}+6241t^{182}7742t^{181}2401t^{180}+1,`$
$`Ans1:=u(t)={\displaystyle \frac{85}{t^{68}}}+79t^{91}+49t^{90};`$
$`DE2:={\displaystyle \frac{du(t)}{dt}}=`$ $`u^2(t)pm\mathrm{cos}(t)\mathrm{sin}^{(m1)}(t)+p^2\mathrm{sin}^{2m}(t),`$
$`Ans2:=u(t)=p\mathrm{sin}^m(t);`$
$`DE3:={\displaystyle \frac{du(t)}{dt}}=`$ $`u^2(t)+pt^{m1}me^{nt}+pt^mne^{nt}+ql{\displaystyle \frac{\mathrm{ln}^{l1}(t)}{t}}`$
$`t^{2m}p^2e^{2nt}2t^mpe^{nt}q\mathrm{ln}^l(t)q^2\mathrm{ln}^{2l}(t),`$
$`Ans3:=u(t)=pt^me^{nt}+q\mathrm{ln}^l(t).`$
Let us make some remarks here. For specific equation we can transform operator solution to improve, e.g. rate of convergence. In the cases under consideration in this section it is very useful with help of operator identities (see below (10) (12)) to modify the operator solution to the following form:

$`u(t)=`$ $`[{\displaystyle \frac{}{\omega }}\mathrm{exp}\{{\displaystyle }dt\omega [{\displaystyle \frac{da(t)}{dt}}a^2(t)]\}\widehat{T}\mathrm{exp}\{{\displaystyle }dt\omega \times `$
$`\times \mathrm{exp}\{{\displaystyle }dt\omega [{\displaystyle \frac{da(t)}{dt}}a^2(t)]\times `$
$`\times {\displaystyle \frac{^2}{\omega ^2}}\mathrm{exp}\{{\displaystyle }dt\omega [{\displaystyle \frac{da(t)}{dt}}a^2(t)]\}e^{\omega C}]_{\omega =0}.`$
Experimentation with prototype of the algorithm leads us to conclusion that the  *algorithm is workable* and can solve (at least in sense of particular solutions) a set of equations which Maple cannot solve now. It is really the *solving* algorithm as against from algorithms based on certain solvability conditions or mapping into separable equations.

Central failure of this algorithm consists in the fact that it *heavily relies on basic CAS procedures and requires noticeable resources*.

It is clear that we can try to find *some* particular solutions (under different values of $`c`$) or even *general* solution of DE by the same way and exploit any form of operator solutions. However the summation of functional series with several variables is for existent CAS substantially more resources consumimg task.

In next sections we examine some ways for solution transformations, which would allow to simplify procedure of solutions calculation.

## 5 Basic transformation identities

Our aim is to reduce the equation solution from operator form (3) to more practical, calculable or usual (standard) expression. To do so it is necessary to find ways to evaluate chronological operator exponents or convert it to combinations of some operators with relatively simple properties, which would permit to carry out all the operators and obtain as a result the solution of equation ”by quadratures”.

One of the main tool kits in the operator method is a family of identities for chronological operator exponents for solutions transformations.

At present time the most-used method to ”disentangle” of operator exponents is well-known Campbell-Baker-Hausdorff (CBH) expansion, which for chronological exponents has the following version:

$`\widehat{T}_0`$ $`\mathrm{exp}\{{\displaystyle _a^\tau }d\xi \widehat{A}(\xi )\}\widehat{B}(\tau )\widehat{T}\mathrm{exp}\{{\displaystyle _a^\tau }d\xi \widehat{A}(\xi )=\widehat{B}(\tau )+{\displaystyle _a^\tau }d\tau _1[\widehat{B}(\tau ),\widehat{A}(\tau _1)]+\mathrm{}`$
$`\mathrm{}+{\displaystyle _a^\tau }𝑑\tau _1{\displaystyle _a^{\tau _1}}𝑑\tau _2\mathrm{}{\displaystyle _a^{\tau _{n1}}}𝑑\tau _n[[\mathrm{}[[\widehat{B}(\tau ),\widehat{A}(\tau _1)],\widehat{A}(\tau _2)]\mathrm{}],\widehat{A}(\tau _n)]+`$
$`+{\displaystyle _a^\tau }d\tau _1{\displaystyle _a^{\tau _1}}d\tau _2\mathrm{}{\displaystyle _a^{\tau _n}}d\tau _{n+1}\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^{\tau _{n+1}}}d\xi \widehat{A}(\xi )\}\times `$
$`\times [[\mathrm{}[[\widehat{B}(\tau ),\widehat{A}(\tau _1)],\widehat{A}(\tau _2)]\mathrm{}],\widehat{A}(\tau _{n+1})]\widehat{T}\mathrm{exp}\{{\displaystyle _a^{\tau _{n+1}}}𝑑\xi \widehat{A}(\xi )\},`$
where $`[\widehat{B}(\tau ),\widehat{A}(\tau _1)]=\widehat{B}(\tau )\widehat{A}(\tau _1)\widehat{A}(\tau _1)\widehat{B}(\tau )`$ is the standard notation of operators commutator.

For purpose of this paper the following related identities are of fundamental importance (all operators are supposed to be linear):

$`\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau \widehat{A}(\tau )\}\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau \widehat{B}(\tau )\}=`$
$`=\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau [\widehat{A}(\tau )+\widehat{T}\mathrm{exp}\{{\displaystyle _a^\tau }𝑑\xi \widehat{A}(\xi )\}\widehat{B}(\tau )\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^\tau }𝑑\xi \widehat{A}(\xi )\}]\}`$ (10)
or in more general form

$`\widehat{A}(t)`$ $`\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau \widehat{B}(\tau )\}=`$
$`=\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}d\tau [({\displaystyle \frac{\widehat{A}(\tau )}{\tau }})\widehat{A}^1(\tau )+\widehat{A}(\tau )\widehat{B}(\tau )\widehat{A}^1(\tau )\}\widehat{A}(a)`$ (11)
end

$`\widehat{T}\mathrm{exp}\{`$ $`{\displaystyle _a^t}d\tau [\widehat{A}(\tau )+\widehat{B}(\tau )]\}=\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}d\tau \widehat{A}(\tau )\}\times `$
$`\times \widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau \widehat{T}_0\mathrm{exp}\{{\displaystyle _a^\tau }𝑑\xi \widehat{A}(\xi )\}\widehat{B}(\tau )\widehat{T}\mathrm{exp}\{{\displaystyle _a^\tau }𝑑\xi \widehat{A}(\xi )\}\}.`$ (12)
Correctness of these identities is easily verified by differentiations (the left and right sides satisfy the same differential equations with the same initial conditions).

From identity (12) follows the identity for chronological exponents with reverse ordering:

$`\widehat{T}_0\mathrm{exp}\{`$ $`{\displaystyle _a^t}d\tau [\widehat{A}(\tau )+\widehat{B}(\tau )]\}=`$
$`=\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}d\tau \widehat{T}_0\mathrm{exp}\{{\displaystyle _a^\tau }d\xi \widehat{A}(\xi )\}\widehat{B}(\tau )\widehat{T}\mathrm{exp}\{{\displaystyle _a^\tau }d\xi \widehat{A}(\xi )\}\}\times `$
$`\times \widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau \widehat{A}(\tau )\}.`$ (13)
We will need later an variant of identity of type (11), when we do not demand from stationary operator $`\widehat{D}`$ to be reversible:

$`\widehat{D}\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau \widehat{B}(\tau )\widehat{D}\}=\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau \widehat{D}\widehat{B}(\tau )\}\widehat{D}.`$ (14)
The identity (10) imparts group properties to chronological operators. We do not evolve here this interesting matter.

## 6 Generalized shift operators

To transform formal solutions to ordinary expressions we have to express explicitly the action of corresponding linear operators. Every solved equation gives us an example of definite action of given linear operator. So the list of ”good” operators is nonempty. Our goal here is to study the properties of certain class of operators that we will name as generalized shift operators.

The Taylor expansion for sufficiently arbitrary function $`\mathrm{\Phi }(x)`$ can be expressed as

$$\mathrm{\Phi }(x+\alpha )=\underset{k=0}{\overset{\mathrm{}}{}}\frac{1}{k!}\alpha ^k\frac{d^k}{dx^k}\mathrm{\Phi }(x)=\mathrm{exp}\{\alpha \frac{d}{dx}\}\mathrm{\Phi }(x).$$

Here $`\alpha `$ can be either constant or function of arguments, which are different from $`x`$ (in the second case $`\frac{d}{dx}`$ is replaced by $`\frac{}{x}`$ . If we read this expression from the right side, we can obtain the result of action of exponential form of *shift operator* $`\mathrm{exp}\{\alpha \frac{d}{dx}\}`$ on function $`\mathrm{\Phi }(x)`$ :

$$\mathrm{exp}\{\alpha \frac{d}{dx}\}\mathrm{\Phi }(x)=\mathrm{\Phi }(\mathrm{exp}\{\alpha \frac{d}{dx}\}x)=\mathrm{\Phi }(x+\alpha )$$

or

$$\mathrm{exp}\{\alpha (t)\frac{}{x}\}\mathrm{\Phi }(x)==\mathrm{\Phi }(\mathrm{exp}\{\alpha (t)\frac{}{x}\}x)=\mathrm{\Phi }(x+\alpha (t)).$$
(15)
Let us consider shift operators with new variable $`y=\psi (x)`$ . It is obvious, that

$$\mathrm{exp}\{\alpha (t)\frac{}{y}\}\mathrm{\Phi }(y)=\mathrm{exp}\{\alpha (t)\frac{}{\psi (x)}\}\mathrm{\Phi }(\psi (x))=$$

$$=\mathrm{exp}\{\alpha (t)\frac{1}{\psi (x)^{}}\frac{}{x}\}\mathrm{\Phi }(\psi (x))=\mathrm{\Phi }(\psi (x)+\alpha (t)).$$

When $`\psi (x)`$ is fixed function, the set of shift operators with different parameters $`\alpha (t)`$ forms Abelian group. Shift operators with different $`\psi (x)`$ do not commute.

It is quite obvious that all examined shift operators correspond to DEs with separable variables.

Let us consider now one important property of the operators of the following type

$$\widehat{T}\mathrm{exp}\{_a^t𝑑\tau f(\tau ,x)\frac{}{x}\},$$

which we will call as generalized shift operators.

Let $`g(t,x)`$ and $`f(t,x)`$ are ordinary functions (linear operators of multiplication by value of function). Since the commutator

$$[g(t,x),f(\tau ,x)\frac{}{x}]=f(\tau ,x)(\frac{g(t,x)}{x})$$

is ordinary function, so from CBH - expansion we can conclude that following expression

$$G(t,x)=\widehat{T}\mathrm{exp}\{_a^t𝑑\tau f(\tau ,x)\frac{}{x}\}g(t,x)\widehat{T}_0\mathrm{exp}\{_a^t𝑑\tau f(\tau ,x)\frac{}{x}\}$$

is ordinary function too and we can compute it simply as

$$G(t,x)=\widehat{T}\mathrm{exp}\{_a^t𝑑\tau f(\tau ,x)\frac{}{x}\}g(t,x).$$
(16)
We find that

$`\widehat{T}`$ $`\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau f(\tau ,x){\displaystyle \frac{}{x}}\}g_1(t,x)g_2(t,x)=`$
$`=\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}d\tau f(\tau ,x){\displaystyle \frac{}{x}}\}g_1(t,x)\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}d\tau f(\tau ,x){\displaystyle \frac{}{x}}\}\times `$
$`\times \widehat{T}\mathrm{exp}\{{\displaystyle _a^t}d\tau f(\tau ,x){\displaystyle \frac{}{x}}\}g_2(t,x)=`$
$`=[\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau f(\tau ,x){\displaystyle \frac{}{x}}\}g_1(t,x)][\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau f(\tau ,x){\displaystyle \frac{}{x}}\}g_2(t,x)]=`$
$`=G_1(t,x)G_2(t,x).`$ (17)
If now expand $`g(t,x)`$ into the Taylor power series with respect to $`x`$, use (17) and then collect series backwards, we obtain analogous to (15) very important property (homomorphy):

$$G(t,x)=\widehat{T}\mathrm{exp}\{_a^t𝑑\tau f(\tau ,x)\frac{}{x}\}g(t,x)=g(t,\widehat{T}\mathrm{exp}\{_a^t𝑑\tau f(\tau ,x)\frac{}{x}\}x).$$
(18)
Now we are ready to consider another variants of operator solutions for equation (5) and farther modifications of solving algorithm.

## 7 Solution for first-order ODE in form of generalized shift operator

Let us start from operator solution of equation (5) in form (7) but let us rewrite it for convenience in following notation

$$S(t,c,\omega )=\widehat{T}\mathrm{exp}\{_a^t𝑑\tau \omega f(\tau ,\frac{}{\omega })\}e^{\omega c},$$

where $`c=u(t)|_{t=a}`$. With help of identity (11) we find

$$S(t,c,\omega )=e^{\omega c}\widehat{T}\mathrm{exp}\{_a^t𝑑\tau \omega e^{\omega c}f(\tau ,\frac{}{\omega })e^{\omega c}\}1.$$

CHB-expantion gives

$$e^{\omega c}f(\tau ,\frac{}{\omega })e^{\omega c}=f(\tau ,\frac{}{\omega })+c[f(\tau ,\frac{}{\omega }),\omega ]+\frac{c^2}{2}[[f(\tau ,\frac{}{\omega }),\omega ],\omega ]+\mathrm{}$$
(19)
Since it is easily proved by mathematical induction that

$$[\frac{^m}{\omega ^m},\omega ]=m\frac{^{m1}}{\omega ^{m1}},$$

then from expantion of $`f(x)`$ into power series and backwards summation it follows that

$$[f(\frac{}{\omega }),\omega ]=f^{}(\frac{}{\omega })(f^{}(x)=\frac{f(x)}{x}),$$

therefore we can conclude with taking into account (19) and property of the shift operator (15) that

$$e^{\omega c}f(\tau ,\frac{}{\omega })e^{\omega c}=f(\tau ,\frac{}{\omega })+cf^{}(\tau ,\frac{}{\omega })+\frac{c^2}{2}f^{\prime \prime }(\tau ,\frac{}{\omega })+\mathrm{}=$$

$$=f(\tau ,c+\frac{}{\omega })=\mathrm{exp}\{\frac{}{\omega }\frac{}{c}\}f(\tau ,c)\mathrm{exp}\{\frac{}{\omega }\frac{}{c}\}.$$

Then

$`S(t,c,\omega )`$ $`=e^{\omega c}\mathrm{exp}\{{\displaystyle \frac{}{\omega }}{\displaystyle \frac{}{c}}\}\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau (\omega {\displaystyle \frac{}{c}})f(\tau ,c)\}\mathrm{exp}\{{\displaystyle \frac{}{\omega }}{\displaystyle \frac{}{c}}\}1=`$
$`=e^{\omega c}\mathrm{exp}\{{\displaystyle \frac{}{\omega }}{\displaystyle \frac{}{c}}\}\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau [\omega f(\tau ,c){\displaystyle \frac{f(\tau ,c)}{c}}f(\tau ,c){\displaystyle \frac{}{c}}]\}1.`$
By expanding the chronological exponent with help of identity (12) and properties (16), (17) of generalized shift operators we obtain that

$`S(t,c,\omega )`$ $`=e^{\omega c}\mathrm{exp}\{{\displaystyle \frac{}{\omega }}{\displaystyle \frac{}{c}}\}\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}d\tau f(\tau ,c){\displaystyle \frac{}{c}}\}\times `$
$`\times \mathrm{exp}\{{\displaystyle _a^t}𝑑\tau g(\tau ,c)\}\mathrm{exp}\{\omega {\displaystyle _a^t}𝑑\tau G(\tau ,c)\},`$ (20)
where

$$g(\tau ,c)=\widehat{T}_0\mathrm{exp}\{_a^\tau 𝑑\xi f(\xi ,c)\frac{}{c}\}\frac{f(\tau ,c)}{c}$$

and

$$G(\tau ,c)=\widehat{T}_0\mathrm{exp}\{_a^\tau 𝑑\xi f(\xi ,c)\frac{}{c}\}f(\tau ,c).$$

If we differentiate the last expression with respect to $`\omega `$ we find that

$`{\displaystyle \frac{S(t,c,\omega )}{\omega }}`$ $`=cS(t,c,\omega )+e^{\omega c}\mathrm{exp}\{{\displaystyle \frac{}{\omega }}{\displaystyle \frac{}{c}}\}\widehat{T}\mathrm{exp}\{{\displaystyle _a^t}d\tau f(\tau ,c){\displaystyle \frac{}{c}}\}\times `$
$`\times \mathrm{exp}\{{\displaystyle _a^t}𝑑\tau g(\tau ,c)\}\mathrm{exp}\{\omega {\displaystyle _a^t}𝑑\tau G(\tau ,c)\}{\displaystyle _a^t}𝑑\tau G(\tau ,c).`$
If in the second item to recover initial operator form we receive

$$\frac{S(t,c,\omega )}{\omega }=cS(t,c,\omega )+S(t,c,\omega )_a^t𝑑\tau G(\tau ,c),$$

inasmuch as $`S(t,c,\omega )|_{\omega =0}1`$, then

$$u(t,c)=\frac{S(t,c,\omega )}{\omega }|_{\omega =0}=c+_a^t𝑑\tau G(\tau ,c),$$

whence it follows that

$$u(t,c)=c+_a^t𝑑\tau \widehat{T}_0\mathrm{exp}\{_a^\tau 𝑑\xi f(\xi ,c)\frac{}{c}\}f(\tau ,c)$$

or finally

$$u(t,c)=\widehat{T}_0\mathrm{exp}\{_a^t𝑑\tau f(\tau ,c)\frac{}{c}\}c.$$
(21)
All above consideration is nonrigourous enough. So it is very important to verify our conclusion by direct substitution of obtained solution into equation (5). After differentiation of (21) and with property (18) we can be sure that (21) is really the solution of the equation (5).

So we have obtained solution of first-order equation in more simple operator form. We can prepare now a prototype of solving algorithm based on solution (21) by the way similar to considered in Section 4. It would be easily found that in this case algorithm can handle more complicated equations and with more efficiency. But it remains here (in sense of present CAS abilities) the complication in obtaining general solutions too. In following Section we consider possibilities to obtain general solution out of known particular solution.

Let us finish this Section with brief remark about interconnection between first-order ODE and linear first-order PDE.

If we introduce the function

$$z(t,c)=\widehat{T}\mathrm{exp}\{_a^t𝑑\tau f(\tau ,c)\frac{}{c}\}c$$
(22)
it is obvious that it satisfies equation

$$\frac{z(t,c)}{t}+f(t,c)\frac{z(t,c)}{c}=0.$$

If we act now by operator

$$\widehat{T}_0\mathrm{exp}\{_a^t𝑑\tau f(\tau ,c)\frac{}{c}\}$$

on both sides of (22), we get

$$z(t,u(t,c))=c$$

– the well-known equation expressing interconnection between solutions of first-order ODE and linear first-order PDE.

Besides that that we *obtain* (operator) solutions of the equations, a nuance here consists in fact that conventional derivation of this interconnection is carried out by means of geometrical considerations (method of characteristics). We use only algebraic transformations and this approach is applicable for more complicated problems. We have the chance to ascertain it in the next Section.

## 8 General solution in terms of one particular solution

If we are able to find some particular solutions of equation (5), then arises the question: can we express the general solution of the equation through these particular solutions? One of the possible answers is the famous solution of this problem by Lie . In terms of operator method Lie’s result is: the general solution can be expressed by relatively simple function of set of particular solutions iff operator solution (7) can be expressed as a product of simple shift operators of type (15). On the other hand it is well-known that, e.g., if we know only one particular solution for Riccati equation the general solution can be found through quadratures.

Let us analyse this question from operator method point of view.

From solution (21) and with help of identities (11), (14) we find that

$`{\displaystyle \frac{u(t,c)}{c}}`$ $`=\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau {\displaystyle \frac{}{c}}f(\tau ,c)\}1=`$
$`=\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau [{\displaystyle \frac{f(\tau ,c)}{c}}f(\tau ,c)+f(\tau ,c){\displaystyle \frac{}{c}}]\}1=`$
$`=\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau \widehat{T}_0\mathrm{exp}\{{\displaystyle _a^\tau }𝑑\xi f(\xi ,c){\displaystyle \frac{}{c}}\}{\displaystyle \frac{f(\tau ,c)}{c}}\}.`$
If denote

$$\frac{f(\tau ,c)}{c}=\mathrm{\Phi }(\tau ,c),$$

then on account of (21) we can find that $`u(t,c)`$ satisfies the following integro-differential equation:

$$\frac{u(t,c)}{c}=\mathrm{exp}\{_a^t𝑑\tau \mathrm{\Phi }(\tau ,u(\tau ,c))\}.$$
(23)
Solution of such type of equations is the challenge.

Let us introduce the following functional

$$S[c;\omega (\xi )]=\mathrm{exp}\{_{\mathrm{}}^{\mathrm{}}𝑑\xi \omega (\xi )u(\xi ,c)\}.$$

Since

$$\frac{S[c;\omega (\xi )]}{c}=S[c;\omega (\xi )]_{\mathrm{}}^{\mathrm{}}𝑑\xi \omega (\xi )\frac{u(\xi ,c)}{c},$$

then to the integro-differential equation (23) corresponds linear equation for $`S[c,\omega (\xi )]`$ with variational derivatives

$$\frac{S[c;\omega (\xi )]}{c}=_{\mathrm{}}^{\mathrm{}}𝑑\xi \omega (\xi )\mathrm{exp}\{_a^\xi 𝑑\tau \mathrm{\Phi }(\tau ,\frac{\delta }{\delta \omega (\tau )})\}S[c;\omega (\xi )]$$

with initial condition

$$S[c;\omega (\xi )]|_{c=c_0}=\mathrm{exp}\{_{\mathrm{}}^{\mathrm{}}𝑑\xi \omega (\xi )u_0(\xi )\},u_0(t)=u(t,c)|_{c=c_0}.$$

Its operator solution is

$`S[c;\omega (\xi )]=`$
$`=\widehat{T}_c\mathrm{exp}\{{\displaystyle _{c_0}^c}𝑑\sigma {\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\xi \omega (\xi )\mathrm{exp}\{{\displaystyle _a^\xi }𝑑\tau \mathrm{\Phi }(\tau ,{\displaystyle \frac{\delta }{\delta \omega (\tau )}})\}\}\mathrm{exp}\{{\displaystyle _{\mathrm{}}^{\mathrm{}}}𝑑\xi \omega (\xi )u_0(\xi )\}.`$
It is obvious that solution in this form is very intricate. So we have to simplify it. If we repeat the chain of transformations almost literally like that in preceding Section, we obtain the solution in more simple form

$$u(t,c)=\mathrm{exp}\{(cc_0)_{\mathrm{}}^{\mathrm{}}𝑑\xi \mathrm{exp}\{_a^\xi 𝑑\tau \mathrm{\Phi }(\tau ,u_0(\tau ))\}\frac{\delta }{\delta u_0(\xi )}\}u_0(t).$$
(24)
Now we have to examine that (24) is really the solution of equations (23) and (5). The verification is not difficult matter here.

So we conclude that if given first-order ODE has general solution it can be expressed in terms of one particular solution of the equation. Certainly we can find the well-known result for Riccati equation from (24).

The expression (24) represents the algorithm for computing the general solution from known particular solution. We have not tried to experiment with this algorithm yet.

## 9 General solution in terms of arbitrary function and solvability conditions

Let now $`f(t,c)=g(t,c)+h(t,c)`$, then we can find the following form of general solution of equation (5)

$`u(t,c)`$ $`=\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau f(\tau ,c){\displaystyle \frac{}{c}}\}c=\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau [g(\tau ,c)+h(\tau ,c)]{\displaystyle \frac{}{c}}\}c=`$
$`=\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}d\tau [\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^\tau }d\xi h(\xi ,c){\displaystyle \frac{}{c}}\}g(\tau ,c)]\times `$
$`\times \mathrm{exp}\{{\displaystyle _a^\tau }d\xi [\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^\xi }d\zeta h(\zeta ,c){\displaystyle \frac{}{c}}\}{\displaystyle \frac{h(\xi ,c)}{c}}]\}{\displaystyle \frac{}{c}}\}\times `$
$`\times \widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau h(\tau ,c){\displaystyle \frac{}{c}}\}c`$
or

$`u(t,c)`$ $`=\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}d\tau [\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^\tau }d\xi h(\xi ,c){\displaystyle \frac{}{c}}\}[f(\tau ,c)h(\tau ,c)]]\times `$
$`\times \mathrm{exp}\{{\displaystyle _a^\tau }d\xi [\widehat{T}_0\mathrm{exp}\{{\displaystyle _a^\xi }d\zeta h(\zeta ,c){\displaystyle \frac{}{c}}\}{\displaystyle \frac{h(\xi ,c)}{c}}]\}{\displaystyle \frac{}{c}}\}\times `$
$`\times \widehat{T}_0\mathrm{exp}\{{\displaystyle _a^t}𝑑\tau h(\tau ,c){\displaystyle \frac{}{c}}\}c,`$ (25)
where we can consider $`h(t,c)`$ as an arbitrary differentiable function.

If it is known the general solution

$$z(t,c)=\widehat{T}_0\mathrm{exp}\{_a^t𝑑\tau h(\tau ,c)\frac{}{c}\}c,$$

then with help of property (18) we obtain that

$$u(t,c)=\widehat{T}_0\mathrm{exp}\{_a^t𝑑\tau \frac{f(\tau ,z(\tau ,c))\frac{z(\tau ,c)}{\tau }}{\frac{z(\tau ,c)}{c}}\frac{}{c}\}z(t,c).$$
(26)
As far as $`h(t,c)`$ is an arbitrary function, therefore preceding expression is valid for any differentiable function $`z(t,c)`$. The particular solution of equation (5) with initial condition $`u(t,c)|_{t=a}=c`$ is obtained when $`z(t,c)|_{t=a}=c`$.

In spite of obviousness and simplicity the form (26) has important meaning. And first of all it gives some recipes for finding ordinary (standard) solution of the initial equation (5).

In contrast to all above mentioned operator forms of solution the last described is not algorithmical as long as $`z(t,c)`$ remains to be an arbitrary function. We do not analyse completely here all consequences of interpretations of solution in form (26). Let restrict ourselves by consideration of some solvability conditions only.

Having noted that if

$$\mathrm{\Phi }(t,c)=\frac{f(t,z(t,c))\frac{z(t,c)}{t}}{\frac{z(t,c)}{c}}$$
(27)
corresponds to right-hand part of solvable first-order ODE, then the general solution of given equation can be obtained too.

If we require $`\mathrm{\Phi }(t,c)=0`$, we need to solve the same given equation. In another cases we will need to solve *different* equation and have chance really to solve it. Main difficulties on such a way arise in the fact that we can not circumscribe all solved DE by single expression, so we have the dilemma: if we choose a certain $`z(t,c)`$ we are not sure that the equation with corresponding $`\mathrm{\Phi }(t,c)`$ is solvable and vice versa – we are not sure that the equation (27) in respect to $`z(t,c)`$ with some $`\mathrm{\Phi }(t,c)`$ selected in advance is solvable.

The situation becomes simpler if we content ourselves with some *solvability conditions*. It is clear that we are interested to pose the broader solvability conditions.

To formulate an algorithmical solvability conditions we have to fix some finite set of functions $`z(t,c)`$ – the set of trial functions. Then by using some criterion we successively examine these function trying to find at least one of them (generally speaking, there are a family of such functions) that satisfy the criterion. If such a function is found, then solution of given equation is obtained almost automatically by substitution of this function to (26).

The following important classes of equations that permit obtaining of general solutions as a whole can illustrate this approach.

*Linear equations*. Here $`\mathrm{\Phi }(t,c)=\alpha (t)c+\beta (t)`$. It is obvious that solvability condition in this case can be formulated as follows. If $`z(t,c)`$ is found to be satisfy the criterion

$$\frac{^2}{c^2}\{\frac{f(t,z(t,c))\frac{z(t,c)}{t}}{\frac{z(t,c)}{c}}\}=0,$$
(28)
then the given equation can be solved through the function $`z(t,c)`$ and a linear equation. Since in this case $`u(t,c)=A(t)z(t,c)+B(t)`$ then here we use linear transformations of variable $`u(t,c)`$. Such type of solving method is due to Liouville. Relatively recent implementations of this approach to CAS make a good showing in solving many of Kamke examples .

*Equations with separable variables*. Here $`\mathrm{\Phi }(t,c)=\alpha (t)\beta (c)`$. The solvability condition in this case can be formulated as follows. If $`z(t,c)`$ is found to be satisfy the criterion

$$\frac{^2}{tc}\{\mathrm{ln}\{\frac{f(t,z(t,c))\frac{z(t,c)}{t}}{\frac{z(t,c)}{c}}\}\}=0,$$
(29)
then the given equation can be solved through the function $`z(t,c)`$ and an equation with separable variables. The suitable functions $`\alpha (t)`$ and $`\beta (c)`$ under known $`z(t,c)`$ are easily obtained from (27).

Reduction of an equation to the one with separable variables of course is the well-known method, but as a matter of fact for solving non-linear equations it is used mainly heuristically simply by guessing appropriate substitution. In resembling a way L. Chan and E.S. Cheb-Terrab have implemented recently so-colled *hypergeometric solutions* method  where the sets of functions $`z(t,c)`$ and substitutions are fixed and only second-order linear (or Riccati) equations are considered. In contrast to their method we here have to fix *only* the set of trial functions $`z(t,c)`$ so our approach can handle broader family of equations.

*Bernoulli equations*. Here $`\mathrm{\Phi }(t,c)=\alpha (t)c^\eta +\beta (t)c`$ ($`\eta 1`$). The solvability condition in this case can be formulated as follows. If $`z(t,c)`$ is found to be satisfy the criterion

$$\frac{}{c}\{\frac{c\frac{}{c}\{\frac{f(t,z(t,c))\frac{z(t,c)}{t}}{\frac{z(t,c)}{c}}\}\frac{f(t,z(t,c))\frac{z(t,c)}{t}}{\frac{z(t,c)}{c}}}{c^\eta }\}=0,$$
(30)
then the given equation can be solved through the function $`z(t,c)`$ and an Bernoulli equation. The suitable functions $`\alpha (t)`$ and $`\beta (t)`$ under known $`z(t,c)`$ are obtained from (27).

Despite the fact that transformations here are simple enough we do not know any real solving method based on this scheme.

We can easily find examples of equations that are solved in such a way and are not solved say by Maple.

It is obvious that list (28)-(30) is far from to be complete, analogous solvability conditions can be obtained from another solved equations. And what is more the solution (26) can be iterated what lead to schemes with successive variables changing (substitutions). In other words we can use combinations of conditions too. The general strategy for such type of algorithms needs further development.

## 10 Conclusions

We presented some ways in solving DE by operator method. Besides first-order DE we can obtain operator solutions for linear and non-linear equations of arbitrary order. The properties of its solutions expect further studies.

In conclusion we believe that we succeeded in demonstration the fact that operator forms of DE solutions can be handled analytically none the worse than ordinary functions. In some cases its transformation properties are more comfortable than, e.g., for some special functions.