Work done is given by $W = F_ids_i$ and in angular displacement moment form it is $W = \tau_i d\theta_i$.
I want to get the work done in using quaternion instead. Like
$W = \tau_i q_i$
I cannot figure out a way to do a dot product between moment and angular displacement since quaternion have four elements. Of course I can keep it in form of angles and change it to quaternion but I don't want that.
The total work done would be
$$ W = \int (\boldsymbol{\tau} \cdot \boldsymbol{\omega}) {\rm d}t $$
I think you are asking about the case of a free-floating object whose axis of rotation changes with time.
In this case the rotational velocity vector is found from the vector part of the changing quaternion in the following fashion
$$\begin{aligned}\dot{q} & =\tfrac{1}{2}\begin{pmatrix}\boldsymbol{\omega}\\ 0 \end{pmatrix}q\\ \begin{pmatrix}\boldsymbol{\omega}\\ 0 \end{pmatrix} & =2\dot{q}q^{-1}\\ \begin{pmatrix}\boldsymbol{\omega}\\ 0 \end{pmatrix} & =2\begin{pmatrix}\dot{\boldsymbol{q}}_{v}\\ \dot{q}_{s} \end{pmatrix}\begin{pmatrix}-\boldsymbol{q}_{v}\\ q_{s} \end{pmatrix}\\ \begin{pmatrix}\boldsymbol{\omega}\\ 0 \end{pmatrix} & =2\begin{bmatrix}q_{s}+[\boldsymbol{q}_{v}\times] & -\boldsymbol{q}_{v}\\ \boldsymbol{q}_{v}^{\intercal} & q_{S} \end{bmatrix}\begin{pmatrix}\dot{\boldsymbol{q}}_{v}\\ \dot{q}_{s} \end{pmatrix} \\ \begin{pmatrix}\boldsymbol{\omega}\\ 0 \end{pmatrix}&=2\begin{pmatrix}q_{s}\boldsymbol{q}_{v}-\boldsymbol{q}_{v}\dot{q}_{s}+\boldsymbol{q}_{v}\times\dot{\boldsymbol{q}}_{v}\\ q_{s}\dot{q}_{s}+\boldsymbol{q}_{v}\cdot\dot{\boldsymbol{q}}_{v} \end{pmatrix}\\ \end{aligned}$$
Above the $[\boldsymbol{q}_v\times]$ notation represents a 3×3 matrix called the cross product operator.
Note that $q$ is the orientation of the body (quaternion) and $\dot{q}$ is the rate of change of the orientation.
$$\boldsymbol{\omega}=2\left(q_{s}\boldsymbol{q}_{v}-\boldsymbol{q}_{v}\dot{q}_{s}+\boldsymbol{q}_{v}\times\dot{\boldsymbol{q}}_{v}\right)$$
The problem is that you cannot just integrate the power relationship to get to work since the rotational components are not necessarily constant (in general).
You can work through the math to get
$$\begin{aligned}{\rm d}W & =\left(\boldsymbol{\tau}\cdot\boldsymbol{\omega}\right){\rm d}t\\ & =2\left(\boldsymbol{\tau}\cdot\left(q_{s}\dot{\boldsymbol{q}}_{v}-\boldsymbol{q}_{v}\dot{q}_{s}+\boldsymbol{q}_{v}\times\dot{\boldsymbol{q}}_{v}\right)\right){\rm d}t\\ & =2\left(q_{s}\left(\boldsymbol{\tau}\cdot\dot{\boldsymbol{q}}_{v}\right)-\dot{q}_{s}\left(\boldsymbol{\tau}\cdot\boldsymbol{q}_{v}\right)+\boldsymbol{\tau}\cdot\left(\boldsymbol{q}_{v}\times\dot{\boldsymbol{q}}_{v}\right)\right){\rm d}t \end{aligned}$$
and with the shorthand notation
$$\begin{aligned}{\rm d}\boldsymbol{q}_{v} & =\dot{\boldsymbol{q}}_{v}{\rm d}t & {\rm d}q_{s} & =\dot{q}_{s}{\rm d}t\end{aligned}$$
$$ \boxed{ {\rm d}W=2\left(q_{s}\left(\boldsymbol{\tau}\cdot{\rm d}\boldsymbol{q}_{v}\right)-\left(\boldsymbol{\tau}\cdot\boldsymbol{q}_{v}\right){\rm d}q_{s}+\left(\boldsymbol{\tau}\times\boldsymbol{q}_{v}\right)\cdot{\rm d}\boldsymbol{q}_{v}\right) } $$
Note that with a fixed axis of rotation $\boldsymbol{\hat{z}}$, you have
$$ \begin{aligned}q & =\begin{pmatrix}\boldsymbol{\hat{z}}\sin\tfrac{\theta}{2}\\ \cos\tfrac{\theta}{2} \end{pmatrix}\\ \boldsymbol{\omega} & =\boldsymbol{\hat{z}}\dot{\theta}\\ \boldsymbol{\tau} & =\boldsymbol{\hat{z}}\tau \end{aligned}$$
and
$$ {\rm d}W = (\boldsymbol{\tau} \cdot \boldsymbol{\omega}) {\rm d}t = \tau \dot{\theta} {\rm d}t = \tau {\rm d}\theta $$
which is the relationship you initially mentioned.
