The reason behind this is not a mathematical reason but rather is an attempt to give an application to the Laplace transform to the analysis of physical systems. This is because any real physical system must necessarily be a casual system.
Straight to the point, these are the answers to your questions:
- The Laplace transform is used because it is more generic and provide more information than the Fourier transform. Furthermore, it is easier to manipulate.
The use of unilateral or bilateral transform should be done with extreme care, depending on the type of causality of the system being analyzed:
- No-Causal System (Anticipative System): Bilateral definition (complete definition):
$$
\begin{align}
\mathcal{L}\left\{ f(t) \right\} &\triangleq \int_{-\infty}^{+\infty} e^{-st} f(t)dt = F(s)
&&&
\mathcal{F}\left\{ f(t) \right\} &\triangleq \int_{-\infty}^{+\infty} e^{-j\omega t} f(t)dt = F(\omega)
\end{align}
$$
- Causal System (Non-anticipative System): Unilateral definition:
$$
\begin{align}
\mathcal{L}\left\{ f(t) \right\} &\triangleq \int_{0}^{+\infty} e^{-st} f(t)dt = F(s)
&&&
\mathcal{F}\left\{ f(t) \right\} &\triangleq \int_{0}^{+\infty} e^{-j\omega t} f(t)dt = F(\omega)
\end{align}
$$
The choice of using the Fourier transform instead of the Laplace transform, is fully valid. But remember three key things:
- Fourier provides less information than Laplace.
- Fourier is more complex than Laplace.
- Use a bilateral or unilateral Fourier definition, according to the causality of the system.
If you do not know or fully understand what I'm talking about, let me explain...
The Black Box Model
Suppose you have a black box in which a signal $x(t)$ is inputted. This signal is processed in the black box and produces an output signal $y(t)$. Then it says that this black box processes the input signal and it produces the output signal, both in the time domain. Note that the black box model can be used to model any type of system. There are no limitations on it.
This black box can be modeled mathematically with a function $h(t)$ (also in the time domain) so that you can set the output signal $y(t)$ by convolution process of the system function $h(t)$ and input signal $x(t)$. Algebraically:
$$
y(t) = h(t) \ast x(t)
$$
How the Black Box Model is related to the Laplace and Fourier transforms?
The relationship between the Laplace/Fourier transform and convolution is the following property:
$$
y(t) = h(t) \ast x(t) \quad\xrightarrow{\mathcal{F}}\quad Y(\omega) = H(\omega) X(\omega)
$$
Is interpreted as the transform of a convolution between $h(t)$ and $x(t)$ is multiplication of $H(\omega)$ and $X(\omega)$. So, the convolution in the time domain becomes a multiplication in the frequency domain. This is very useful since compute a multiplication is much simpler than computing a convolution.
Note that there is no restriction on using the Fourier transform or the Laplace transform to do this. However, the definition of the Laplace transform is more generic with respect to the Fourier transform, since $s=\sigma + \omega i$.
$$
\begin{align}
\mathcal{L}\left\{ f(t) \right\} &\triangleq \int_{-\infty}^{+\infty} e^{-st} f(t)dt = F(s)
&&&
\mathcal{F}\left\{ f(t) \right\} &\triangleq \int_{-\infty}^{+\infty} e^{-j\omega t} f(t)dt = F(\omega)
\end{align}
$$
$$
\Longrightarrow\quad
\begin{matrix}
f(t) & \quad\xrightarrow{\mathcal{L}}\quad & F(s) \\
f(t) & \quad\xrightarrow{\mathcal{F}}\quad & F(\omega)
\end{matrix}
\quad\therefore\quad \left. F(s) \right|_{\sigma = 0} = F(\omega)
$$
For this reason, the Laplace transform is used. In addition, the variable $ s $ provides additional information because:
- Variable $\omega$ is related to the permanent system response; and
- Variable $\sigma$ is related to the inherent attenuation of the system.
$$
y(t) = h(t) \ast x(t) \quad\xrightarrow{\mathcal{L}}\quad Y(s) = H(s) X(s)
$$
Therefore, the Laplace transform is a powerful tool to analyze a general system, and even allows obtain the function $h(t)$ and/or $y(t)$ using the inverse transform with the following procedure:
How does all this is related with unilateral or bilateral definition of the Laplace transform?
There is a property of physical systems called causality, where the output of the black box depends on past and current input but not future inputs. This idea is intuitive. A clear example is the following: the weather tomorrow will depend of present and past conditions, but can never depend on the weather the day after tomorrow. In the Wikipedia link, this idea is better explained.
And what influences a system is causal or not? That it influences the definition of the system function $h(t)$:
$$
h(t) \mbox{ is a casual system} \quad\Leftrightarrow\quad \forall t\in\mathbb{R},\, t < 0:\quad h(t) = 0
$$
Therefore, when computing the Laplace transform of a causal system, you get the unilateral definition of the Laplace transform:
$$
\mathcal{L}\left\{ h(t) \right\} \triangleq \underbrace{\int_{-\infty}^{+\infty} e^{-st} h(t)dt}_{\mbox{Bilateral Def.}} = \underbrace{\int_{0}^{+\infty}e^{-st} h(t)dt }_{\mbox{Unilateral Def.}}
$$
Note that it is extremely important force analysis of a physical system to a causal system, especially since the system function $h(t)$ can be mathematically defined for $t <0$. For this reason, it applied directly the unilateral Laplace transform for real physical systems (or more generally speaking, for any casual system).
There are also no-causal systems (such as software that digitally processes an image) where the system state may depend on future states of the system ("future" pixels). In these cases, you must necessarily apply the bilateral transform.
This answer assumes the following definitions of the Fourier, Laplace, and two-sided Laplace transforms.
$$\mathcal{F}_t[f(t)](\omega)=\int\limits_{-\infty}^{\infty } f(t)\,e^{-i \omega t}\,dt\tag{1}$$
$$\mathcal{L}_t[f(t)](s)=\int\limits_0^{\infty} f(t)\,e^{-s t}\,dt\tag{2}$$
$$\mathcal{BL}_t[f(t)](s)=\int\limits_{-\infty}^{\infty} f(t)\,e^{-s t}\,\tag{3}dt$$
Note the Fourier transform defined in formula (1) above is equivalent to the bilateral Laplace transform defined in formula (3) above evaluated at $s=i \omega$.
The Fourier and Laplace transforms of $\cos(\omega_0\,t)$ are as follows:
$$\mathcal{F}_t[\cos(\omega_0\,t)](\omega)=\int\limits_{-\infty}^{\infty}\cos (\omega_0\,t)\,e^{-i \omega t}\,dt=\pi\,\delta(\omega+\omega_0)+\pi\,\delta(\omega-\omega_0),\quad\omega\in\mathbb{R}\tag{4}$$
$$\mathcal{L}_t[\cos(\omega_0\,t)](s)=\int\limits_0^{\infty} \cos(\omega_0\,t)\,e^{-s t}\,dt=\frac{s}{s^2+\omega_0^2},\quad\Re(s)>0\tag{5}$$
The Fourier and bilateral Laplace transforms of $\cos(\omega_0\,t)\,\theta(t)$ are
$$\mathcal{F}_t[\cos(\omega_0\,t)\,\theta(t)](\omega)=\int\limits_{-\infty}^{\infty}\cos(\omega_0\,t)\,\theta(t)\,e^{-i \omega t}\,dt$$ $$=\int_0^{\infty}\cos(\omega_0\,t)\,e^{-i \omega t}\,dt=\frac{\pi}{2} \delta(\omega+\omega_0)+\frac{\pi}{2} \delta(\omega-\omega_0)+\frac{i \omega}{\omega_0^2-\omega^2},\quad\omega\in\mathbb{R}\tag{6}$$
$$\mathcal{BL}_t[\cos(\omega_0\,t)\,\theta(t)](s)=\int\limits_{-\infty}^{\infty}\cos(\omega_0\,t)\,\theta(t)\,e^{-s t}\,dt=\int\limits_0^{\infty} \cos(\omega_0\,t)\,e^{-s t}\,dt=\frac{s}{s^2+\omega_0^2},\ \Re(s)>0\tag{7}$$
where
$$\theta(t)=\left\{
\begin{array}{cc}
0 & t<0 \\
1 & t>0 \\
\end{array}\right.\tag{8}$$
is the Heaviside step function.
Consider the following limit representation of the bilateral Laplace transform of $\cos(\omega_0\,t)$.
$$\mathcal{BL}_t[\cos(\omega_0\,t)](s)=\underset{T\to\infty}{\text{lim}}\left(\int\limits_{-T}^{T}\cos(\omega_0\,t)\,e^{-s t}\,dt\right)=\underset{T\to\infty}{\text{lim}}\left(\frac{2 (s \sinh(s T) \cos(T \omega_0)+\omega_0\cosh (s T) \sin (T \omega_0))}{s^2+\omega_0^2}\right)\tag{9}$$
Substituting $T=2 \pi f$ and $s=i \omega$ in formula (9) above and accounting for the resulting equivalence to formula (4) above leads to the limit representation
$$\delta(\omega+\omega_0)+\delta(\omega-\omega_0)=\underset{f\to\infty}{\text{lim}}\left(\frac{2}{\pi} \frac{\omega \cos(2 \pi f \omega_0) \sin (2 \pi f \omega )-\omega_0 \sin(2 \pi f \omega_0) \cos(2 \pi f \omega )}{\omega ^2-\omega_0^2}\right),\quad\omega\in\mathbb{R}\tag{10}$$
which is equivalent to the limit representation
$$\delta(\omega+\omega_0)+\delta(\omega-\omega_0)=\underset{f\to\infty}{\text{lim}}\left(2 f\,\text{sinc}(2 \pi f (\omega+\omega_0))+2 f\,\text{sinc}(2 \pi f (\omega -\omega_0))\right),\quad\omega\in\mathbb{R}.\tag{11}$$
The relevant question seems to be whether $\mathcal{F}_t[\cos(\omega_0\,t)\,\theta(t)](\omega)$ can be approximated by $\underset{\epsilon\to 0^+}{\text{lim}}\left(\mathcal{BL}_t[\cos(\omega_0\,t)\,\theta(t)](\epsilon+i \omega)\right)=\underset{\epsilon\to 0^+}{\text{lim}}\left(\mathcal{L}_t[\cos(\omega_0\,t)](\epsilon+i \omega)\right)$ for $\omega\in\mathbb{R}$ which, based on formulas (5) to (7) above, is equivalent to the following.
$$\frac{\pi}{2} \delta(\omega+\omega_0)+\frac{\pi}{2} \delta(\omega-\omega_0)+\frac{i \omega}{\omega_0^2-\omega^2}\approx\underset{\epsilon\to 0^+}{\text{lim}}\left(\frac{\epsilon+i \omega}{(\epsilon+i \omega)^2+\omega_0^2}\right),\quad\omega\in\mathbb{R}\tag{12}$$
Figure (1) below illustrates the imaginary part of the right side of formula (12) in orange overlaid on the imaginary part of the left-side of formula (12) in blue where both sides are evaluated at $\omega_0=1$ and the right side is evaluated at $\epsilon=10^{-6}$. Note for $\epsilon=0$ the right-side of formula (12) above becomes equivalent to the last term on the left-side of formula (12) above (which is the imaginary part of the left-side of formula (12) above), so clearly $\mathcal{F}_t[\cos(\omega_0\,t)\,\theta(t)](\omega)\ne\mathcal{L}_t[\cos(\omega_0\,t)](i \omega)$.
Figure (1): Illustration of imaginary parts of right and left sides of formula (12) (orange overlaid on blue) evaluated at $\omega_0=1$ and $\epsilon=10^{-6}$
Figures (2) below illustrates the real part of the right-side of formula (12) above evaluated at $\omega_0=1$ and $\epsilon=10^{-6}$. Note the evaluation seems to suggest the real part of the right-side of formula (12) above may perhaps approximate a limit representation for the real-part of the left side of formula (12) above (i.e. the two Dirac delta terms).
Figure (2): Illustration of real part of right-side of formula (12) evaluated at $\omega_0=1$ and $\epsilon=10^{-6}$
The analysis and results illustrated above seem to suggest the left-side of formula (12) above can be approximated by the limit representation on the right-side of formula (12) above, but the following integral which assumes $\omega_0=1$ provides further evidence.
$$\int\limits_{-2}^2 \left(\frac{\pi}{2} \delta(\omega+1)+\frac{\pi}{2} \delta(\omega-1)+\frac{i \omega}{1-\omega^2}\right)\,d\omega =\pi\approx \underset{\epsilon\to 0^+}{\text{lim}}\left(\frac{i}{2} (\log(\epsilon (\epsilon-4 i)-3)-\log(\epsilon (\epsilon+4 i)-3))\right)=\underset{\epsilon\to 0^+}{\text{lim}}\left(\int\limits_{-2}^2 \frac{\epsilon+i \omega}{(\epsilon+i \omega)^2+1}\,d\omega\right)\tag{13}$$
Note for the integral of the imaginary part of the left-side of formula (12) above the Cauchy principal value $\mathcal{P}\int\limits_{-2}^2 \frac{i \omega}{1-\omega^2}\,d\omega=0$ is consistent with the fact that $\frac{i \omega}{1-\omega^2}$ is an odd function of $\omega$ as illustrated in Figure (1) above.
For $\epsilon=10^{-6}$ the right side of formula (13) above evaluates to $3.14159\,+0.i$ which is very close to the value of $\pi$.
Best Answer
What the question means, I think, is to take the Fourier Transform in the spatial variable and the Laplace Transform in the time variable.
First, taking the Laplace Transform $\mathcal{L}[u(x,t](s) = w(x, s)$ and letting $u_0(x) := u(x, 0)$ be a known initial condition, we obtain:
$$sw(x,s) - u_0(x) -\kappa\frac{\partial^2 w(x,t)}{\partial x^2} = S_0 \delta(x)$$
which is equivalent to
$$\left(s - \kappa\frac{\partial^2}{\partial x^2}\right)w(x,s) = S_0\delta(x) + u_0(x)$$
Taking the Fourier Transform of this equation, letting $\mathcal{F}[v(x, s)](k) = v(k, s)$, we obtain
$$\left(s+\kappa k^2\right)v(k, s) = S_0 + \int_{-\infty}^{\infty}u_0(x)e^{-ikx}dx$$
or, equivalently
$$v(k,s) = \frac{1}{s + \kappa k^2}\left(S_0 + \int_{-\infty}^{\infty}u_0(x)e^{-ikx}dx\right)$$
(depending on your Fourier Transform convention, various factors of ${2\pi}^{-1}$ may appear).
Now, we only have to invert the two transforms, and then we are done.
The spatial stuff looks pretty messy, but the Laplace Transform can be done easily, as we only have one $s$ in there.
In fact, we have a function that looks something like this:
$$v(k, s) = \frac{F(k)}{s + a(k)}$$
which has an Inverse Laplace Transform of $p(k, t) = \mathcal{L}^{-1}[v(k,s)](t) = F(k)\exp(-a(k)t)$, so we are left with
$$p(k, t) = \exp\left(-\kappa k^2 t\right)\left(S_0 + \int_{-\infty}^{\infty}u_0(x)e^{-ikx}dx\right)$$
The Inverse Fourier Transform of the first term can be reduced to calculating a Gaussian integral, giving us the typical Gaussian solution one would expect from the Heat Equation as a homogeneous solution. This calculation is outlined below:
\begin{align} u_H(x, t) &= \mathcal{F}[S_0 \exp\left(-\kappa k^2 t\right)](x, t)\\ &= \frac{S_0}{2\pi}\int_{-\infty}^\infty \exp(-\kappa t k^2 + i x k)dk\\ &= \frac{S_0}{2\pi}\int_{-\infty}^\infty \exp\left(-\kappa t\left(k^2 - \frac{ixk}{\kappa t}\right)\right)dk\\ &= \frac{S_0}{2\pi}\int_{-\infty}^\infty\exp\left(-\kappa t\left(k-\frac{ix}{2\kappa t}\right)^2 - \frac{x^2}{4\kappa t}\right)dk\\ &= \frac{S_0}{2\pi}\exp\left(-\frac{x^2}{4\kappa t}\right)\int_{-\infty}^\infty \exp\left(-\kappa t\left(k-\frac{ix}{2\kappa t}\right)^2\right)dk \end{align}
Letting $z := k - \frac{ix}{2\kappa t}$, we get that $dk = dz$. The contour on integration shifts to a line parallel to the real axis, but it is not too difficult to show that if one were to consider a a rectangular closed contour the runs along both this line and the real axis, the contribution of the short edges of the rectangle would go to zero as the length of the line approaches infinity. Thus, we can shift the integral back onto the real axis, as the integral along the real axis and along the line parallel to it must be equal in value.
\begin{align} u_H(x, t) &= \frac{S_0}{2\pi}\exp\left(-\frac{x^2}{4\kappa t}\right)\int_{-\infty}^\infty \exp\left(-\kappa tz^2\right)dz\\ &= \frac{S_0}{2\pi}\exp\left(-\frac{x^2}{4\kappa t}\right)\sqrt{\frac{\pi}{\kappa t}}\\ &= \frac{S_0}{\sqrt{4\pi\kappa t}}\exp\left(-\frac{x^2}{4\kappa t}\right) \end{align}
This is the known solution to the Heat Equation as a homogeneous solution.
The inhomogeneous solution can only be written down in integral form, as we do not know anything about the initial distribution $u_0(x)$. It is given by:
\begin{align} u_I(x, t) &= \frac{1}{2\pi}\int_{-\infty}^\infty\int_{-\infty}^{\infty}\exp\left(-\kappa k^2 t\right)e^{ikx} e^{-iky} u_0(y) dy dk\\ &= \frac{1}{2\pi}\int_{-\infty}^\infty u_0(y)\int_{-\infty}^\infty \exp\left(-\kappa t k^2 + ik(x-y)\right)dkdy\\ &= \frac{1}{\sqrt{4\pi \kappa t}}\int_{-\infty}^\infty u_0(y)\exp\left(-\frac{(x-y)^2}{4\kappa t}\right) dy \end{align}
which is preciously the convolution of the initial condition $u_0(x)$ with the fundamental homogeneous solution $u_H(x, t)$.
The final solution is then given as the sum of the homogeneous and the inhomogeneous solution.
EDIT: The above was done for a general $u_0(x)$, because I missed that $u_0(x)$ had been specified in the questions. Given that we have $u_0(x) = \delta(x)$, we can calculate the final integral easily:
\begin{align} u_I(x,t) &= \frac{1}{\sqrt{4\pi \kappa t}}\int_{-\infty}^\infty \delta(y)\exp\left(-\frac{(x-y)^2}{4\kappa t}\right) dy\\ &= \frac{1}{\sqrt{4\pi \kappa t}}\exp\left(-\frac{x^2}{4\kappa t}\right) \end{align}
This gives us the full solution