Solution 1:

The problem with the initial approach is that along $C_2$, $\dfrac{e^{iaz}}{x^3+1}$ does not vanish as $R\to\infty$. For $a\gt0$, it vanishes in the upper half-plane, but blows up exponentially in the lower half-plane.


Convert from Oscillating to Monotonic Integral

Let's compute the integral of $\dfrac{e^{iaz}}{z^3+1}$ over $\gamma$, where $$ \gamma=[0,R]\cup Re^{i[0,\pi/2]}\cup[iR,0]\tag{1} $$ The integral of $\dfrac{e^{iaz}}{z^3+1}$ along $Re^{i[0,\pi/2]}$ is less than $\frac{R\pi/2}{R^3-1}$ which vanishes as $R\to\infty$.

The residue of $\dfrac{e^{iaz}}{z^3+1}$ at $z=e^{i\pi/3}$ is $\frac13e^{-a\sqrt3/2}\left[\cos\left(\frac a2-\frac{2\pi}3\right)+i\sin\left(\frac a2-\frac{2\pi}3\right)\right]$.

Since $z=e^{i\pi/3}=\frac12+i\frac{\sqrt3}2$ is the only singularity inside $\gamma$, we get $$ \begin{align} &\frac{2\pi i}3e^{-a\sqrt3/2}\left[\cos\left(\frac a2-\frac{2\pi}3\right)+i\sin\left(\frac a2-\frac{2\pi}3\right)\right]\\ &=\int_\gamma\frac{e^{iaz}}{z^3+1}\mathrm{d}z\\ &=\int_0^\infty\frac{e^{iax}}{x^3+1}\,\mathrm{d}x-i\int_0^\infty\frac{e^{-ax}}{1-ix^3}\,\mathrm{d}x\\ &=\int_0^\infty\frac{e^{iax}}{x^3+1}\,\mathrm{d}x-i\int_0^\infty\frac{e^{-ax}(1+ix^3)}{1+x^6}\,\mathrm{d}x\tag{2} \end{align} $$ Solving for the integral we want: $$ \begin{align} \int_0^\infty\frac{e^{iax}}{x^3+1}\,\mathrm{d}x &=-\int_0^\infty\frac{x^3e^{-ax}}{1+x^6}\,\mathrm{d}x-\frac{2\pi}3e^{-a\sqrt3/2}\sin\left(\frac a2-\frac{2\pi}3\right)\\ &+i\left[\int_0^\infty\frac{e^{-ax}}{1+x^6}\,\mathrm{d}x+\frac{2\pi}3e^{-a\sqrt3/2}\cos\left(\frac a2-\frac{2\pi}3\right)\right]\tag{3} \end{align} $$ Therefore, separating real and imaginary parts, we have $$ \int_0^\infty\frac{\cos(ax)}{x^3+1}\mathrm{d}x =-\int_0^\infty\frac{x^3e^{-ax}}{1+x^6}\,\mathrm{d}x-\frac{2\pi}3e^{-a\sqrt3/2}\sin\left(\frac a2-\frac{2\pi}3\right)\tag{4} $$ and $$ \int_0^\infty\frac{\sin(ax)}{x^3+1}\mathrm{d}x =\int_0^\infty\frac{e^{-ax}}{1+x^6}\,\mathrm{d}x+\frac{2\pi}3e^{-a\sqrt3/2}\cos\left(\frac a2-\frac{2\pi}3\right)\tag{5} $$


Asymptotic Expansion

Equations $(4)$ and $(5)$ show that as $a\to\infty$, the integrals vanish. In fact, it is not difficult to obtain an asymptotic expansion from the right hand side of these integrals: $$ \int_0^\infty\frac{x^3e^{-ax}}{1+x^6}\,\mathrm{d}x \sim\frac{3!}{a^4}-\frac{9!}{a^{10}}+\frac{15!}{a^{16}}-\dots\tag{6} $$ and $$ \int_0^\infty\frac{e^{-ax}}{1+x^6}\,\mathrm{d}x \sim\frac{0!}{a^1}-\frac{6!}{a^7}+\frac{12!}{a^{13}}-\dots\tag{7} $$ Unfortunately, neither Mathematica nor I have found any simpler formulas for these exponential integrals. However, Mathematica does verify that equations $(4)$ and $(5)$ hold numerically for several values of $a$.


Convergent Series

Although equations $(4)$ and $(5)$ are valid for any $a\ge0$ and make the integral easier to compute numerically, and $(6)$ and $(7)$ give good approximations for large $a$, the question has now been amended to express interest in small $a$.

Integration by parts gives $$ \newcommand{\Ein}{\operatorname{Ein}} \begin{align} \int_a^\infty\frac{e^{-x}}{x}\mathrm{d}x &=\int_0^a\frac{1-e^{-x}}{x}\mathrm{d}x-\log(a)-\gamma\\[6pt] &=\Ein(a)-\log(a)-\gamma\tag{8} \end{align} $$ where $$ \begin{align} \Ein(x)&=\sum_{k=1}^\infty(-1)^{k-1}\frac{x^k}{k\,k!}\tag{9}\\ \end{align} $$ or to remove oscillation for positive $x$: $$ \begin{align} e^x\Ein(x)&=\sum_{k=1}^\infty\frac{H_k}{k!}x^k\tag{10}\\ \end{align} $$ Therefore, $$ \begin{align} \int_0^\infty\frac{e^{-x}\mathrm{d}x}{x+a} &=e^a\int_a^\infty\frac{e^{-x}\mathrm{d}x}{x}\\[6pt] &=e^a\left(\Ein(a)-\log(a)-\gamma\right)\tag{11} \end{align} $$ A change of variables, partial fractions yields, and applications of $(10)$ and $(11)$ yield $$ \begin{align} &\int_0^\infty\frac{x^3e^{-ax}}{1+x^6}\,\mathrm{d}x\\ &=a^2\int_0^\infty\frac{x^3e^{-x}}{a^6+x^6}\,\mathrm{d}x\\ &=\small-\frac16\int_0^\infty\left(\frac1{x+ia}+\frac\omega{x+i\omega a}+\frac{\omega^2}{x+i\omega^2a}+\frac1{x-ia}+\frac\omega{x-i\omega a}+\frac{\omega^2}{x-i\omega^2a}\right)e^{-x}\,\mathrm{d}x\\ &=-\frac13\int_0^\infty\mathrm{Re}\left(\frac1{x+ia}+\frac\omega{x+i\omega a}+\frac{\omega^2}{x+i\omega^2a}\right)e^{-x}\,\mathrm{d}x\\ &=\sum_{k=0}^\infty(-1)^k\frac{H_{6k+2}}{(6k+2)!}a^{6k+2}\\ &+\small\frac13\left[\cos(a)-\cos\left(\tfrac12a\right)\cosh\left(\tfrac{\sqrt3}2a\right)-\sqrt3\sin\left(\tfrac12a\right)\sinh\left(\tfrac{\sqrt3}2a\right)\right](\log(a)+\gamma)\\ &+\frac\pi{18}\sin\left(\tfrac12a\right)\left[3\cosh\left(\tfrac{\sqrt3}2a\right)-2\sinh\left(\tfrac{\sqrt3}2a\right)\right]\\ &+\frac\pi{6\sqrt3}\cos\left(\tfrac12a\right)\left[2\cosh\left(\tfrac{\sqrt3}2a\right)-3\sinh\left(\tfrac{\sqrt3}2a\right)\right]\\ &-\frac\pi6\sin(a)\tag{12} \end{align} $$ and $(12)$ provides a real method to compute the integral appearing in $(4)$.

Similarly for the integral appearing in $(5)$ $$ \begin{align} &\int_0^\infty\frac{e^{-ax}}{1+x^6}\,\mathrm{d}x\\ &=a^5\int_0^\infty\frac{e^{-x}}{a^6+x^6}\,\mathrm{d}x\\ &=\small\frac16\int_0^\infty\left(\frac{i}{x+ia}+\frac{i\omega}{x+i\omega a}+\frac{i\omega^2}{x+i\omega^2a}-\frac{i}{x-ia}-\frac{i\omega}{x-i\omega a}-\frac{i\omega^2}{x-i\omega^2a}\right)e^{-x}\,\mathrm{d}x\\ &=-\frac13\int_0^\infty\mathrm{Im}\left(\frac1{x+ia}+\frac\omega{x+i\omega a}+\frac{\omega^2}{x+i\omega^2a}\right)e^{-x}\,\mathrm{d}x\\ &=\sum_{k=0}^\infty(-1)^{k+1}\frac{H_{6k+5}}{(6k+5)!}a^{6k+5}\\ &+\small\frac13\left[\sin(a)+\sin\left(\tfrac12a\right)\cosh\left(\tfrac{\sqrt3}2a\right)-\sqrt3\cos\left(\tfrac12a\right)\sinh\left(\tfrac{\sqrt3}2a\right)\right](\log(a)+\gamma)\\ &+\frac\pi{18}\cos\left(\tfrac12a\right)\left[3\cosh\left(\tfrac{\sqrt3}2a\right)-2\sinh\left(\tfrac{\sqrt3}2a\right)\right]\\ &-\frac\pi{6\sqrt3}\sin\left(\tfrac12a\right)\left[2\cosh\left(\tfrac{\sqrt3}2a\right)-3\sinh\left(\tfrac{\sqrt3}2a\right)\right]\\ &+\frac\pi6\cos(a)\tag{13} \end{align} $$

Solution 2:

$\newcommand{\angles}[1]{\left\langle\, #1 \,\right\rangle} \newcommand{\braces}[1]{\left\lbrace\, #1 \,\right\rbrace} \newcommand{\bracks}[1]{\left\lbrack\, #1 \,\right\rbrack} \newcommand{\ceil}[1]{\,\left\lceil\, #1 \,\right\rceil\,} \newcommand{\dd}{{\rm d}} \newcommand{\ds}[1]{\displaystyle{#1}} \newcommand{\dsc}[1]{\displaystyle{\color{red}{#1}}} \newcommand{\expo}[1]{\,{\rm e}^{#1}\,} \newcommand{\fermi}{\,{\rm f}} \newcommand{\floor}[1]{\,\left\lfloor #1 \right\rfloor\,} \newcommand{\half}{{1 \over 2}} \newcommand{\ic}{{\rm i}} \newcommand{\iff}{\Longleftrightarrow} \newcommand{\imp}{\Longrightarrow} \newcommand{\Li}[1]{\,{\rm Li}_{#1}} \newcommand{\pars}[1]{\left(\, #1 \,\right)} \newcommand{\partiald}[3][]{\frac{\partial^{#1} #2}{\partial #3^{#1}}} \newcommand{\pp}{{\cal P}} \newcommand{\root}[2][]{\,\sqrt[#1]{\vphantom{\large A}\,#2\,}\,} \newcommand{\sech}{\,{\rm sech}} \newcommand{\sgn}{\,{\rm sgn}} \newcommand{\totald}[3][]{\frac{{\rm d}^{#1} #2}{{\rm d} #3^{#1}}} \newcommand{\ul}[1]{\underline{#1}} \newcommand{\verts}[1]{\left\vert\, #1 \,\right\vert}$ \begin{align}&\color{#66f}{\large% \int_{0}^{\infty}{\cos\pars{ax} \over x^{3} + 1}\,\dd x} =a^{2}\int_{0}^{\infty}{\cos\pars{x} \over x^{3} + \verts{a}^{3}}\,\dd x =-\,{1 \over 3\verts{a}^{3}}\sum_{n\ =\ -1}^{1}x_{n}\int_{0}^{\infty} {\cos\pars{x} \over x - x_{n}}\,\dd x \\[5mm]&\mbox{where}\qquad x_{-1} \equiv \exp\pars{-\,{\pi\ic \over 3}}\verts{a}\,,\quad x_{0} \equiv -\verts{a}\,,\quad x_{1} \equiv \exp\pars{{\pi\ic \over 3}}\verts{a} \end{align}

Then, \begin{align}&\color{#66f}{\large% \int_{0}^{\infty}{\cos\pars{ax} \over x^{3} + 1}\,\dd x} ={1 \over 3a^{2}}\dsc{\int_{0}^{\infty}{\cos\pars{x} \over x + \verts{a}}\,\dd x} -{2 \over 3\verts{a}^{3}}\,\Re\bracks{x_{1}\dsc{\int_{0}^{\infty} {\cos\pars{x} \over x - x_{1}}\,\dd x}}\quad\pars{1} \end{align}

However, \begin{align}&\dsc{\int_{0}^{\infty}{\cos\pars{x} \over x + \mu}\,\dd x} =\int_{\mu}^{\infty}{\cos\pars{x - \mu} \over x}\,\dd x \\[5mm]&=\cos\pars{\mu}\int_{\mu}^{\infty}{\cos\pars{x} \over x}\,\dd x +\sin\pars{\mu}\int_{\mu}^{\infty}{\sin\pars{x} \over x}\,\dd x \\[5mm]&=-\cos\pars{\mu}\,{\rm Ci}\pars{\mu} +\sin\pars{\mu}\bracks{{\pi \over 2} - \,{\rm Si}\pars{\mu}} \\[5mm]&=\dsc{{\pi \over 2}\,\sin\pars{\mu} -\cos\pars{\mu}\,{\rm Ci}\pars{\mu} -\sin\pars{\mu}\,{\rm Si}\pars{\mu}}\tag{2} \end{align} where $\ds{\,{\rm Ci}}$ and $\ds{\,{\rm Si}}$ are the Cosine Integral and the Sine Integral , respectively.

In using the result $\pars{2}$, the integral $\pars{1}$ can be expressed in terms of the Cosine Integrals and Sine Integrals.