Asymptotic expansion of the Volterra integral equation with series ansatz Problem
I have a problem which can be boiled down to the Volterra integral equation
$$
\begin{aligned}
w(\eta) &\sim \eta \int _0^\infty K(s)F(\eta s)ds
\end{aligned}
$$
for $\alpha \in (\frac 12 , 1)$, with kernel and non-linearity defined as
$$
K(s) := \frac{1}{\Gamma(\alpha)}(1-s)^{\alpha-1} \theta(1-s) \qquad F(\eta s) := (\eta s + \eta_0)^{-1-\alpha} (c_0 + c_1w(\eta s) + c_2 w(\eta s)^2)
$$
where $\theta(s)$ is Heaviside step function is $\eta_0>0$ is a positive constant.
It is assumed that for our choice of coefficients $c_0, c_1, c_2$ the function $w(\eta)\rightarrow \infty$ as $\eta \rightarrow \infty$.
The equation cannot be solved explicitly, the goal is to obtain the asymptotic expansion of $w(\eta)$ as $\eta \rightarrow \infty$.

My attempt
Recall the definition of the Mellin transform
$$
\mathrm M [f(s); z] = \int _0^\infty s^{z-1} f(s) ds.
$$
Once formulated as above, we can use Parseval formula for the Mellin transform and rewrite the problem as
$$
w(\eta) \sim \frac{\eta}{2 \pi i} \int_{c-i \infty}^{c+i \infty} \mathrm{M}[K(s) ; 1-z] \mathrm{M}[F(\eta s) ; z] d z
$$
whereas $\Re (c)$ is in the analiticity strip of both Mellin transforms.
Ansatz
We now make an ansatz
$$
w(\eta) = \sum_{k=1}^{-\infty} d_k \eta^{\alpha k} = d_1 \eta^\alpha + d_0 + d_{-1} \eta ^ {-\alpha} + \dots = d_1 \eta^\alpha + d_0 + O(\eta^{-\alpha})
$$
The idea is to calculate Mellin transforms in the above integral and match the coefficients to calculate $d_1$ and $d_0$.
Mellin transform of the kernel
$$
M[K(s); 1-z] = \frac{\Gamma(1-z)}{\Gamma(1+\alpha-z)} = - \sum _{n=0} ^\infty \frac{(-1)^n}{\Gamma(\alpha-n) n!} \left(\frac{1}{z-(n+1)}\right).
$$
It has poles at $z = 1, 2, 3, \dots$ and hence analytical for $\Re (z) < 1$.
Mellin transform of the non-linearity
If we plug in $w(\eta)$ from ansatz above into the non-linearity $F(\eta s)$, we obtain
$$
\begin{aligned}
F(\eta s) &= (\eta s + \eta _0)^{-1-\alpha}\left\{c_2 d_1^2  (\eta s)^{2\alpha} + (2 c_2 d_0d_1 + c_1 d_1)(\eta s)^{\alpha}\right\} + O((\eta s)^{-1-\alpha} ) \\
&\sim c_2 d_1^2  (\eta s)^{-1-\alpha} + (2 c_2 d_0d_1 + c_1 d_1)(\eta s)^{-1} + O((\eta s)^{-1-\alpha} )
\end{aligned}
$$
Mellin transforms of the first two terms are
$$
\begin{aligned}
M\left[c_2 d_{1}^{2} (\eta s)^{-1+\alpha}; z\right] &= (c_2 d_{1}^{2}) \frac{\eta ^{-z}}{z-(1-\alpha)}\\
M\left[(2 c_2 d_0d_1 + c_1 d_1) (\eta s)^{-1}; z\right] &= (2 c_2 d_0d_1 + c_1 d_1) \frac{\eta ^{-z}}{z-1}.
\end{aligned}
$$
In order for the Mellin transforms to exist, we require $z < 1-\alpha$. Hence $\mathrm{M}[F(\eta s) ; z]$ is analytic for $\Re(z) < 1-\alpha$. So we see that each term in the expansion of $F(\eta s)$ leads to the pole of $\mathrm{M}[F(\eta s) ; z]$ at $z = 1-\alpha, 1, 1+\alpha, \dots$.
Shifting the contour to the right
Recall that using Parseval formula we have rewritten our problem as
$$
w(\eta) \sim \frac{\eta}{2 \pi i} \int_{c-i \infty}^{c+i \infty} \mathrm{M}[K(s) ; 1-z] \mathrm{M}[F(\eta s) ; z] d z.
$$
Now that we obtained analiticity strips of the Mellin transforms above, we see it is required that $\Re (c) < 1 - \alpha$. The idea now is to shift the integration contour to the right and use Cauchy integral formula
$$
f^{(n)}(a)=\frac{n !}{2 \pi i} \oint_{\gamma} \frac{f(z)}{(z-a)^{n+1}} d z
$$
to calculate the integrals around the poles.
Assume at this point that we can shift the contour to the right. We plug in the Mellin transforms calculated above and shift the contour to the right so that $u \in (1, 1+\alpha)$ to obtain first two terms in the integral expansion:
$$
\begin{aligned}
d_1 \eta ^\alpha + d_0 + \dots &\sim \frac{\eta}{2 \pi i} \int_{\gamma_1} \left( \frac{\Gamma(1-z)}{\Gamma(1+\alpha-z)}\right) \left(c_{2} d_{1}^{2}\right)  \frac{\eta^{-z}}{z-(1-\alpha)} \\
&+ \frac{\eta}{2 \pi i} \int_{\gamma_0} \left( \frac{\Gamma(1-z)}{\Gamma(1+\alpha-z)}\right) \left(2 c_{2} d_{0} d_{1}+c_{1} d_{1}\right) \frac{\eta^{-z}}{z-1} \\
&+ \frac{\eta}{2 \pi i} \int_{u - i\infty} ^{u + i \infty} \mathbf{M}[K(s) ; 1-z] \mathbf{M}[F(\eta s) ; z] d z
\end{aligned}
$$
such that $\gamma_1$ encloses $z=1-\alpha$ and $\gamma_0$ encloses $z=1$.
First coefficient $d_1$
The contour integral around $z = 1-\alpha$ can be calculated with Cauchy integral formula 
$$
\frac{\eta}{2 \pi i} \int_{\gamma_1} \left(c_{2} d_{1}^{2}\right) \left( \frac{\Gamma(1-z)}{\Gamma(1+\alpha-z)}\right) \frac{\eta^{-z}}{z-(1-\alpha)} = c_2 d_1^2 \frac{\Gamma(\alpha)}{\Gamma(2\alpha)} \eta ^\alpha
$$
Matching the first coefficient leads to 
$$
d_1 \eta ^\alpha = c_2 d_1^2 \frac{\Gamma(\alpha)}{\Gamma(2\alpha)} \eta ^\alpha
\qquad \text{and thus} \qquad
d_1 = \frac{\Gamma(2\alpha)}{c_2 \Gamma(\alpha)}
$$
Second coefficient $d_0$
The problem that arises now is that both kernel and non-linearity have pole at $z=1$. Hence the integrand has the pole of order 2 and Cauchy integral formula will involve derivative of $\eta ^{-z}$, which will yield the term $\log (\eta)$:
$$
\begin{array}{l}
\frac{\eta}{2 \pi i} \int_{\gamma_0} \frac{\Gamma(1-z)}{\Gamma(1+\alpha-z)} M[F(\eta s) ; z] d z \\
= \eta \cdot \text{Res}(\Gamma(1-z), 1) \cdot \text{Res}(M[F(\eta s); z], 1) \cdot \left. [f'(z)] \right| _{z = 1, f(z) = \frac{\eta^{-z}}{\Gamma(1+\alpha-z)}} \\
= (2c_2 d_{0} d_{1}+c_1 d_{1}) \frac{-\log (\eta) \Gamma(\alpha) + \Gamma'(\alpha)}{\Gamma(\alpha)^2}.
\end{array}
$$
Matching the coefficient $d_0$ leads to
$$
d_0 = (2c_2 d_{0} d_{1}+c_1 d_{1}) \frac{\Gamma'(\alpha)}{\Gamma(\alpha)^2}.
$$
Since we know the $d_1$, the above equation can be solved for $d_0$.
However, we are left with the residual term $- \frac{(2c_2 d_{0} d_{1}+c_1 d_{1}) \Gamma(\alpha)}{\Gamma(\alpha)^2} \log (\eta)$ of order $ \log (\eta)$ of order $ \log (\eta)$, which cannot be matched with anything on the left-hand side, since our ansatz did not have a logarithmic term. Hence this ansatz does not work directly and probably should be modified.

Remark
Extending the ansatz to involve the logarithmic term
$$
w(\eta) = d_1 \eta^\alpha + d_\text{log} \log (\eta) + d_0 + d_{-1} \eta ^ {-\alpha} + \dots
$$
will lead to the pole of order 3 in the integrand, and Cauchy integral formula will involve second derivative and lead to the term of order $\log ^2 (\eta)$ on the right-hand side (including $\log ^2 (\eta)$ will lead to $\log ^3 (\eta)$ on the right-hand side and so on).

In conclusion, I wanted to ask whether there is a way to improve this ansatz so that all coefficients can be matched. On the other hand, maybe some alternative approach to the problem can be considered.
Any help appreciated, thank you.
 A: Definitions
Define $P(x)=c_2 x^2+c_1 x+c_0$.
Rewrite the problem as
$$w(\eta)\sim\frac1{\Gamma(\alpha)}\int^\eta_0 ds\frac{(1-s/\eta)^{\alpha-1}}{(s+\eta_0)^{\alpha+1}}P\circ w(s)$$
As a shorthand, define the operator $$I=\int^\eta_0 ds\frac{(1-s/\eta)^{\alpha-1}}{(s+\eta_0)^{\alpha+1}}$$
Also define $A$ as the 'asymptotic operator', such that $A[f]$ returns the asymptotics of $f$ up to $o(1)$, according to certain asymptotic scale.
Then, the problem becomes a fixed-point problem: $$w=\frac1{\Gamma(\alpha)}\cdot A\circ I\circ P\circ w\qquad(\star)$$

Properties of $I$
It is assumed that $w(\eta)$ grows at most polynomially as $\eta\to\infty$.
First, we investigate $A\circ I$'s action on a monomial $x^q$ ($q\ge0$).
$$\begin{align}
I[s^q](\eta)
&=\int^\eta_0\frac{(1-s/\eta)^{\alpha-1}}{(s+\eta_0)^{\alpha+1}}s^q ds \\
&=\int^\eta_0\frac{s^q}{(s+\eta_0)^{\alpha+1}}ds+
\underbrace{\int^\eta_0\frac{(1-s/\eta)^{\alpha-1}-1}{(s+\eta_0)^{\alpha+1}}s^q ds}_{I'(q)} \\
\end{align}
$$
For $I'(q)$,
$$\begin{align}
I'(q)
&=\int^{\eta}_0 \frac{(1-s/\eta)^{\alpha-1}-1}{(s+\eta_0)^{\alpha+1}}s^q ds \\
&=\sum^\infty_{r=1}\binom{\alpha-1}{r}(-1)^r \eta^{-r}\int^\eta_0\frac{s^{q+r}}{(s+\eta_0)^{\alpha+1}}ds \\
&=\sum^\infty_{r=1}\binom{\alpha-1}{r}(-1)^r \eta^{-r}\left(\frac{\eta^{q+r-\alpha}}{q+r-\alpha}+O(\eta^{q+r-\alpha-1})\right) \\
&=\eta^{q-\alpha}\sum^\infty_{r=1}\binom{\alpha-1}{r}\frac{(-1)^r}{q+r-\alpha}+O(\eta^{q-\alpha-1}) \\
&=\eta^{q-\alpha}\left[\frac{\Gamma(\alpha)\Gamma(q-\alpha)}{\Gamma(q)}-\frac1{q-\alpha}\right]+O(\eta^{q-\alpha-1}) \\
\end{align}
$$
Note that $\gamma(q):=\frac{\Gamma(\alpha)\Gamma(q-\alpha)}{\Gamma(q)}-\frac1{q-\alpha}$ converges for $\Re (q)\ge 0$.

Summary $$I[s^q](\eta)=\int^\eta_0\frac{s^q}{(s+\eta_0)^{\alpha+1}}ds+\gamma(q)\eta^{q-\alpha}+O(\eta^{q-\alpha-1})$$

By differentiating under the integral sign, it can be further proved that
$$I[s^q \ln^N s](\eta)=
\begin{cases}
O(\eta^{q-\alpha}\ln^N \eta), &q>\alpha \\
O\left(\ln^{N+1}\eta\right), &q=\alpha \\
O(1), &\alpha> q\ge 0
\end{cases}
$$
(Asymptotic expansions up to $o(1)$ are given in the appendix below.)

A different ansatz
It is obvious that, to match the leading order on both sides of $(\star)$, $w(\eta)\sim C\eta^{\alpha}$ where $C=\frac{\Gamma(2\alpha)}{c_2 \Gamma(\alpha)}$, as the OP obtained.
Now, let's make a slightly different ansatz:
$$w(\eta)=C\eta^{\alpha}+\sum^N_{r=0}B_r\ln^r\eta$$
Utilizing asymptotic formulae in the appendix, ($D:=\frac{2\Gamma(2\alpha)}{\Gamma^2(\alpha)}$)
$$\begin{align}
\frac1{\Gamma(\alpha)}\cdot A\circ I\circ P\circ w
&=C\eta^{\alpha}+\frac{DB_N}{N+1}\ln^{N+1}\eta \\
&+D\sum^N_{r=2}\left[\frac{B_{r-1}}{r}+\sum^N_{m=r}\binom{m}{r}\gamma^{(m-r)}(\alpha)B_m\right]\ln^r\eta \\
&+D\left[\frac{c_1\Gamma(\alpha)}{2c_2}+2B_0+\sum^N_{m=1}m\gamma^{(m-r)}(\alpha)B_m\right]\ln\eta \\
&+\text{a messy bunch of constants}
\end{align}
$$
(The messy bunch of constants can be found using the asymptotic formulae in appendix.)
The mismatch of order, i.e. the presence of $\ln^{N+1}\eta$, is expected. However, in the following calculations, we will give evidence hinting at rapid decay of $B_r$, such that

*

*The coefficient of the mismatching order, $\frac{DB_N}{N+1}$, vanishes as $N\to\infty$.

*$w(\eta)=C\eta^{\alpha}+\sum^\infty_{r=0}B_r\ln^r\eta$ converges.

If these two conditions are satisfied, an ansatz with $N=\infty$ is a solution to $(\star)$.

Evidence suggesting $N\to\infty$ solves the problem

Claim $$ B_{N-i}=c_i\cdot\frac{N!}{(N-i)!}B_N \qquad (\blacksquare)$$ where $c_i$ is independent of $N$.

By comparing coefficients,
$$B_r=\frac{B_{r-1}}{r}+\sum^N_{m=r}\binom{m}{r}\gamma^{m-r}(\alpha)B_m$$
$$\implies \frac{B_{r-1}}{r}+(\gamma(\alpha)-1)B_r+\sum^N_{m=r+1}\binom{m}{r}\gamma^{m-r}(\alpha)B_m=0 \qquad (1)$$
Let $r=N-i$. After a re-indexing of the sum,
$$\frac{B_{N-i-1}}{N-i}+(\gamma(\alpha)-1)B_{N-i}+\sum^{i-1}_{m=0}\binom{N-m}{N-i}\gamma^{i-m}(\alpha)B_{N-m}=0$$
Substitute in $B_{N-i}=c_i\cdot\frac{N!}{(N-i)!}B_N$,
$$\begin{align}
&c_{i+1}\frac{N! B_N}{(N-i)!}+(\gamma(\alpha)-1)c_i\frac{N! B_N}{(N-i)!} \\
&+\sum^{i-1}_{m=0}\frac{(N-m)!}{(N-i)!(i-m)!}\cdot \gamma^{i-m}(\alpha)\cdot c_m\frac{N!B_N}{(N-m)!} = 0
\end{align}
$$
Here we notice mass cancellations, leading to
$$c_{i+1}=(1-\gamma(\alpha))c_i-\sum^{i-1}_{m=0}\frac{\gamma^{(i-m)}(\alpha)}{(i-m)!}c_m$$
Absence of $N$ supports our claim.
Furthermore, by a contour integration approach, it can be proved that
$$\frac{\gamma^{(x)}(\alpha)}{x!}=\frac{(-1)^{x}}{\Gamma(1-\alpha)}\sum^{\infty}_{n=1}\frac{\Gamma(n+\alpha+1)}{n! n^{x+1}}$$ for $x\ge 1$.
This can be used to estimate $c_i\approx R^i$ for some constant $R>1$.
Thus, we conclude that:

*

*Setting $i\to N-i$ in $(\blacksquare)$, we have $B_i=\frac{c_{N-i}}{i!}N!B_N$, we get $$\frac{B_{i+1}}{B_{i}}\approx \frac{1/R}{i+1}\to 0 \qquad\text{as }i\to\infty$$ As a result, it suggests that $w(\eta)=C\eta^{\alpha}+\sum^\infty_{r=0}B_r\ln^r\eta$ converges.

*Setting $i=N$ in $(\blacksquare)$, we get $B_N=\frac{B_0}{c_N\cdot N!}$ Despite not being explicitly stated, the dependence between $B_0$ and $N$ exists. However,
$B_0$ corresponds to the messy bunch of constants above, and it is hard to analyse it analytically. One way out is altering the problem slightly: suppose $A$ extracts terms of order higher than $O(1)$ (instead of $o(1)$), and add a 'perturbation constant' $\lambda$ to $(\star)$ -
$$\lambda+w=\frac{1}{\Gamma(\alpha)}A\circ I\circ P\circ w \qquad (\star ')$$ Then, the right hand side of $(\star')$ has no $O(1)$ term, so $B_0=-\lambda$ and the 'coefficent of mismatching order' becomes approximately $-\frac{D\lambda}{R^N (N+1)!}\to 0$ as $N\to\infty$.

Practically, most Volterra integral equations are solved numerically. Therefore, it is feasible to directly plug in the infinte series and solve for the coefficients numerically. Noteworthily, viewing $B_r$ as unknowns, $(1)$ creates a system of  linear equations, with an almost triangular coefficient matrix and the constant matrix  $\begin{pmatrix} -B_0 & 0 & \cdots & 0 \end{pmatrix}^T$. Solving the system seems like an easy task for computer.
(This explains why $B_0$ (or $\lambda$) cannot be zero: otherwise the system of equations becomes homogeneous, and since the coefficient matrix is invertible, all $B_r=0$. The solution becomes trivial.)
A very special solution exists when $\displaystyle{\eta_0=\left[\frac{c_1^2(1-c_2)-c_0}{c_1\Gamma(\alpha+1)}\right]^{1/\alpha}}$: $$w(\eta)=C\eta^{\alpha}-c_1$$
One approach that I haven't tried is substituting $\eta=e^z$ - the infinte series then becomes a Maclaurin series, whose coefficients may be found (probably) by Cauchy's integral formula or simply by differentiation. I may add on this later.

Maclaurin series approach
Define $v(z)=w(e^z)$.
Substitute in $\eta=e^z$,
$$v(z)=\frac{e^z}{\Gamma(\alpha)}\int^1_0 \frac{(1-s)^{\alpha-1}}{(e^z s+\eta_0)^{\alpha+1}}P\circ v(z+\ln s)ds$$
The ansatz is $$v(z)=\sum^\infty_{i=0}a_i z^i$$
Let's start from the right hand side.
$$\begin{align}
v(z+\ln s)
&=\sum^\infty_{i=0}\sum^i_{j=0}a_i\binom{i}{j}(\ln s)^{i-j}z^j \\
&=\sum^\infty_{j=0}\left[\sum^\infty_{i=j}a_i\binom{i}{j}(\ln s)^{i-j}\right]z^j \\
&=\sum^\infty_{j=0}\sum^\infty_{i=0}a_{i+j}\binom{i+j}{j}(\ln s)^{i}z^j \\
&=\sum^\infty_{j=0}\sum^\infty_{i=0}A[i,j](\ln s)^{i}z^j \\
\end{align}
$$
where $A[x,y]=a_{x+y}\binom{x+y}{y}$.
Viewing this double series as a 'two-dimensional' Taylor series, Cauchy product immediately gives
$$v^2(z+\ln s)=\sum^\infty_{i=0}\sum^\infty_{j=0}(A*A)[i,j](\ln s)^{i}z^j \\$$
where $(A*A)[x,y]=\sum^x_{m=0}\sum^y_{n=0}A[m,n]A[x-m,y-n]$ is the discrete convolution.
Hence,
$$P\circ v(z+\ln s)=c_0+\sum^\infty_{i=0}\sum^\infty_{j=0}\underbrace{\bigg(c_2(A*A)[i,j]+c_1A[i,j]\bigg)}_{B_{ij}}(\ln s)^i z^j$$
Operating $\frac1{\Gamma(\alpha)}\cdot A\circ I$ on $P\circ v(z+\ln s)$, we can finally convert $(\star)$ into
$$\Gamma(\alpha)\sum^\infty_{j=0}a_j z^j
=\frac{c_0\eta_0^{-\alpha}}{\alpha}
+\sum^\infty_{j=0}\sum^\infty_{i=0}B_{ij}\sum^i_{k=0}\binom{i}{k}\int^\infty_0\frac{\ln^{i-k}s}{(s+\eta_0)^{\alpha+1}}ds \cdot z^{j+k}$$
which is suitable for coefficient comparison and solving by computer. At last, $w(\eta)=v(\ln\eta)$.

Appendix
For $\alpha<q\le 2\alpha$,

$$I[s^q \ln^N s](\eta)=\int^\eta_0\frac{s^q\ln^N s}{(s+\eta_0)^{\alpha+1}}ds+\eta^{q-\alpha}\sum^N_{r=0}\binom{N}{r}\gamma^{(m-r)}(q)\ln^r\eta+O(\eta^{q-\alpha-1}\ln^N \eta)$$

For $q=\alpha$,

$$\begin{align}
I[s^\alpha \ln^N s](\eta)
&=\frac{\ln^{N+1}\eta}{N+1}+\eta_0\cdot\frac{\alpha+1}{N+1}\int^\infty_0\frac{s^\alpha\ln^{N+1}s}{(s+\eta_0)^{\alpha+2}}ds \\
&+\sum^N_{r=0}\binom{N}{r}\gamma^{(m-r)}(\alpha)\ln^r\eta+O(\eta^{-1}\ln^N \eta)
\end{align}
$$

For $0\le q<\alpha$,

$$I[s^q \ln^N s](\eta)=\int^\infty_0\frac{s^q\ln^N s}{(s+\eta_0)^{\alpha+1}}ds+O(\eta^{q-\alpha}\ln^N\eta)$$

