In order to familiarize myself with perturbation methods, I've been trying to derive the Lorentz transformations, given by
\begin{align*} x \rightarrow \frac{x + vt}{\sqrt{1 - v^2}} & = (x + vt)(1 + \frac{1}{2}v^2 + \frac{3}{8}v^4 + ...) \\ & = (x + vt + \frac{1}{2}v^2x + \frac{1}{2}v^3t + \frac{3}{8}v^4x + \frac{3}{8}v^5t + ...) \end{align*} \begin{align*} t \rightarrow \frac{t + vx}{\sqrt{1 - v^2}} & = (t + vx)(1 + \frac{1}{2}v^2 + \frac{3}{8}v^4 + ...) \\ & = (t + vx + \frac{1}{2}v^2t + \frac{1}{2}v^3x + \frac{3}{8}v^4t + \frac{3}{8}v^5x + ...) \end{align*}
using a perturbative approach which entails starting with the Gallilean transformations (i.e. $x \rightarrow x + vt$, $t \rightarrow t$) and in each step adding variations $\delta x$ and $\delta t$, both of order $n$ in $v$ and linear in $x$ and $t$, such that the spacetime interval $t^2 - x^2$ is invariant with respect to the transformations up to order $n+1$ in $v$. However, I keep running into the problem that, in each step of the process, some of the coefficients in my expansions for the transformations are undetermined. For example, for $\delta x$ and $\delta t$ of order $2$ in $v$,
$$x \rightarrow x + vt + \delta x$$
$$t \rightarrow t + vx + \delta t$$
$$\delta x = v^2(a_1x + b_1t) + \mathcal{O}(v^3)$$
$$\delta t = v^2(a_2x + b_2t) + \mathcal{O}(v^3)$$
\begin{align*} t^2 - x^2 & = (t + vx + \delta t)^2 - (x + vt + \delta x)^2 \\ & = (t + vx)^2 + 2(t + vx)\delta t - (x + vt)^2 - 2(x + vt)\delta x + \mathcal{O}(v^4) \\ & = t^2 + 2vxt + v^2x^2 + 2v^2(t + vx)(a_2x + b_2t) \\ & \qquad - x^2 - 2vxt - v^2t^2 - 2v^2(x + vt)(a_1x + b_1t) + \mathcal{O}(v^4) \\ & = t^2 + v^2x^2 + 2v^2a_2xt + 2v^2b_2t^2 \\ & \qquad - x^2 - v^2t^2 - 2v^2a_1x^2 - 2v^2b_1xt + \mathcal{O}(v^3) \\ & = t^2 - x^2 + v^2[(2b_2-1)t^2 + 2(a_2-b_1)xt + (1-2a_1)x^2] + \mathcal{O}(v^3) \end{align*}
I end up with $a_1 = b_2 = \frac{1}{2}$ and $a_2 = b_1 \equiv c_1$, which is so far an undetermined coefficient.
I know from the Taylor expansions of the full transformations above that $c_1$ should be equal to $0$. Is there a way to mathematically determine this coefficient?
The coefficients are underdetermined because the appropriate boundary condition has not been applied yet. The $v$ in the expansion is currently only an arbitrary perturbation parameter, and in order for the perturbation expansion to yield the correct expression for the Lorentz transformations, the parameter $v$ must be equated with the physically real boost velocity that it is supposed to represent.
Consider an object stationary at the origin in a reference frame $K$ with coordinates $(t,x)$; the worldline of this object in $K$ is given by the straight line $(t,0)$. Now consider a boost on $K$ by velocity $v$ in the $-x$ direction and denote the new reference frame as $K'$ using coordinates $(t',x')$. The coordinates of $K'$ may be given by a perturbation expansion in $v$ of $(t + vx + \delta t, x + vt + \delta x)$, and the worldline of the object stationary in $K$ is given by $(t + \delta t, vt + \delta x)$ since $x=0$ for all $t$. For two points on the worldline -- given by $P_1$ and $P_2$ in $K$ and $P_1'$ and $P_2'$ in $K'$ respectively -- the slope of the line between the two points measured in a given reference frame gives the velocity of the object according to the given reference frame. In other words,
$$v_K = \frac{\Delta x}{\Delta t} = 0$$ $$v_{K'} = \frac{\Delta x'}{\Delta t'} = \frac{v\Delta t + \delta (\Delta x)}{\Delta t + \delta (\Delta t)}$$
Since, from the perspective of an observer in $K'$, the velocity of an object stationary with respect to $K$ must be $v$ in the $+x$ direction due to the physical principles of relativity, then $v_{K'}$ must be equal to $v$ up to $\mathcal O (v^n)$. For example, up to $\mathcal O (v)$,
$$\require{cancel} v = \frac{v\Delta t + v^2(\frac{1}{2}\cancelto{0}{\Delta x} + c_1\Delta t) + \mathcal O(v^3)}{\Delta t + v^2(c_1\cancelto{0}{\Delta x} + \frac{1}{2}\Delta t) + \mathcal O(v^3)} = \frac{v(\Delta t + vc_1\Delta t + \mathcal O(v^2))}{\Delta t + \mathcal O(v^2)}$$
If the RHS of the previous equation equals $v$ up to $\mathcal O(v)$, then $c_1 = 0$.
This procedure may be used to fully determine all of the coefficients in the expansion and give
$$x \rightarrow \frac{x + vt}{\sqrt{1-v^2}}$$
$$t \rightarrow \frac{t + vx}{\sqrt{1-v^2}}$$
as the appropriate Lorentz transformations.