Just "fill in the blanks." Suppose we use Thevenin for the left-most divider network to find:
$$V_\text{X}=V_\text{IN}\frac{R_2}{R_2+Z_{C_2}}$$ and $$Z_\text{X}=\frac{R_2\,Z_{C_2}}{R_2+Z_{C_2}}$$
With that Thevenin source (whose \$Z_\text{X}\$ is now in series with \$R_1\$), it then follows that:
$$G_s=\frac{V_\text{o}}{V_\text{i}}=\frac{R_2}{R_2+Z_{C_2}}\cdot\frac{Z_{C_1}}{R_1+Z_{C_1}+\frac{R_2\,Z_{C_2}}{R_2+Z_{C_2}}}$$
Of course, \$Z_C=\frac1{s\,C}\$. So, with a little bit of algebra stuff:
$$G_s=\frac{R_2\,C_2\,s}{s^2+\left(\frac1{R_1\,C_1}+\frac1{R_2\,C_2}+\frac1{R_1\,C_2}\right)s+\frac{1}{R_1\,C_1\,R_2\,C_2}}$$
Set \$\alpha=\frac12 \left(\frac1{R_1\,C_1}+\frac1{R_2\,C_2}+\frac1{R_1\,C_2}\right)\$, \$\omega_{_0}=\frac1{\sqrt{R_1\,C_1\,R_2\,C_2}}\$, and create the unitless \$\zeta=\frac{\alpha}{\omega_{_0}}\$. It follows that \$f_{_0}\approx 112.54\:\text{kHz}\$ and \$\zeta\approx 8.132\$.
We want to find \$K\$ such that:
$$G_s=K\:\left[\frac{2\zeta\,\omega_{_0}\,s}{s^2+2\zeta\,\omega_{_0}\,s+\omega_{_0}^2}\right]$$
(If you are curious about the separation of these two factors, examine what happens to the second factor when \$s=j\,\omega_{_0}\$.)
We find here that \$K=\frac{R_2\,C_2}{R_1\,C_1+R_2\,C_2+R_2\,C_1}\$. (You can instead just compute \$K=\sqrt{G_{j\,\omega_{_0}}\:G_{-j\,\omega_{_0}}}\$, but there are less formal methods of getting to the same place.) With your values, I find \$K\approx 0.08695\$ or that the filter peaks at about \$-21.214\:\text{dB}\$.
The denominator is quadratic and the roots are:
$$\begin{align*}\left\{\begin{array}{l}s_1=-\alpha+\sqrt{\alpha^2-\omega_{_0}^2}\\s_2=-\alpha-\sqrt{\alpha^2-\omega_{_0}^2}\end{array}\right.\end{align*}$$
\$\zeta\$ is handy. The following cases arrive (if you look at the square-root term of \$s_1\$ and \$s_2\$ you may note that it can be imaginary or real):
$$\begin{align*}\text{Damping factor conditions}\left\{\begin{array}{l}\zeta = 1 \left(\alpha=\omega_0\right)&&\text{Critically damped}\\\zeta \gt 1 \left(\alpha\gt \omega_0\right)&&\text{Over-damped}\\\zeta \lt 1 \left(\alpha\lt \omega_0\right)&&\text{Under-damped}\\\zeta = 0&&\text{Un-damped}\end{array}\right.\end{align*}$$
In your case, you have a bandpass because \$\zeta\approx 8.132\$ and it must be the over-damped case. So the square-root part of the solution is real and therefore \$s_1\$ and \$s_2\$ are both real (and different from each other.) Here also, the \$s_1\$ and \$s_2\$ poles actually represent your \$\omega_{_\text{L}}\$ and \$\omega_{_\text{H}}\$:
$$\begin{align*}\left\{\begin{array}{l}\omega_{_\text{L}}=-s_1=\omega_{_0}\left(\zeta-\sqrt{\zeta^2-1}\right)=\omega_{_0}\,\zeta\left(1-\sqrt{1-\frac1{\zeta^2}}\right)\\\omega_{_\text{H}}=-s_2=\omega_{_0}\left(\zeta+\sqrt{\zeta^2-1}\right)=\omega_{_0}\,\zeta\left(1+\sqrt{1-\frac1{\zeta^2}}\right)\end{array}\right.\end{align*}$$
Easily computed from the variables earlier developed. (And note that \$\omega_{_\text{L}}\,\omega_{_\text{H}}=\omega_{_0}^2\$.) So \$f_{\text{L}}\approx 6.946\:\text{kHz}\$ and \$f_{\text{H}}\approx 1.82334\:\text{MHz}\$.
Note that the prior transfer function is but one standard way to write it. It's not the only standard form. Another approach is to simply replace \$s\$ with \$j\,\omega\$ (we are assuming that it doesn't spiral out of control or that it doesn't dampen into nothing -- in short, we are assuming \$\sigma=0\$.) Also, since \$s_1\$ and \$s_2\$ are real roots (for a bandpass filter that is over-damped, by definition), we can re-arrange as follows:
$$\begin{align*}
G_s&=-K\:\left[\frac{2\zeta\,\omega_{_0}\,s}{\left(s-s_1\right)\cdot\left(s-s_2\right)}\right]\\\\
&=-K\:\left[\frac{2\zeta\,\omega_{_0}\,j\,\omega}{\left(j\,\omega+\omega_{_\text{L}}\right)\cdot\left(j\,\omega+\omega_{_\text{H}}\right)}\right]\\\\
&=-K\:\left[\frac{2\zeta\,\frac{j\,\omega}{\omega_{_0}}}{\left(1+\frac{j\,\omega}{\omega_{_\text{L}}}\right)\cdot\left(1+\frac{j\,\omega}{\omega_{_\text{H}}}\right)}\right]
\end{align*}$$
The point is that there are different ways to represent the same thing. Your choice will depend on what you want to emphasize. (And, of course, it's worth the time to play a little with the equations to see where they take you.)
(By the way, the above only applies to the case where we are talking about an over-damped bandpass filter. I made some assumptions about the fact that we have real and distinct roots for this last development.)
Your Results
I notice that you used AC 11, not AC 1, in your simulation. You should correct this. The standard way to approach this is to use AC 1 in your simulation's input source.
If you fix this in your simulation, everything else will work out fine. Here are the results from LTspice:
I think that should suffice.
Summary
You have to take into account the loading that one sub-network has upon another. So you need to use the full Thevenin approach that takes into account all of the loading effects.