For the damped harmonic oscillator equation $$\frac{d^2x}{dt^2}+\frac{c}{m}\frac{dx}{dt}+\frac{k}{m}x=0$$ we get that the general solution is $$x(t)=Ae^{-\gamma t}e^{i\omega_d t}+Be^{-\gamma t}e^{-i\omega_d t}$$ where $\gamma = \frac{c}{2m}$ and $ \omega_d=\sqrt{\omega^2-\gamma ^2}$. Using Euler's equation, we can expand this as follows: $$Ae^{-\gamma t}(\cos(\omega _dt)+i\sin(\omega_d t))+Be^{-\gamma t}(\cos(\omega _dt)-i\sin(\omega_d t))$$ $$\Rightarrow e^{-\gamma t}(A+B)\cos(\omega_d t) +e^{-\gamma t}(Ai-Bi)\sin(\omega_d t)$$ But now we are dealing with a physical problem so we only examine the real part which is $e^{-\gamma t}(A+B)\cos(\omega_d t)$. But this does not have any phase difference. Yet textbooks always make the claim that the real part of the solution is $$e^{-\gamma t}(C)\cos(\omega_d t+\phi)$$ where $\phi$ is some arbitrary initial phase. But where does that initial phase come from if the real part of the solution does not have a phase change in it? I understand that $A$ and $B$ themselves need not be real however I do not understand how this fact could ever lead to a non zero initial phase in the real part of the solution.
This issue has bothered me for quite some time now so any help would be immensely appreciated!
Your split into real and imaginary parts is flawed for generic complex coefficients A and B. In general, you may parameterize them as $$ A=a e^{i\theta+i\phi}, \qquad B=b e^{i\theta-i\phi}, $$ for real a,b,θ,φ. It is evident that θ is superfluous, as its corresponding phase factors out of the linear equation, and thus x.
So only φ is meaningful, and it joins up with the dynamical phase $\omega_d t$, providing a new origin for it. Then, with real a,b, the amended real-imaginary split you attempted would be sound, and the conventional result would follow for real C=a+b.