> In [this paper](https://arxiv.org/abs/2103.04967), on page 3, the authors go from the tight binding model w the Peierls substitution
$$ H = \sum_{i,j} \sum_{a,b} t_{a,b}   \exp\left(i \int_{\textbf{R}_{j,b}}^{\textbf{R}_{i,a}} dr'_{\mu} A_{\mu} (\textbf{r}',t ) \right) c_{i,a}^\dagger c_{j,b} $$
to $$ H= H_0 + \sum_{i,j} \sum_{a,b} t_{a,b} (L_\mu^A A_\mu + 1/2 L_{\mu v}^{AA} A_\mu A_v +..) c_{i,a}^\dagger c_{j,b} $$
where $H_0 $ is $H$ without the exponential, and then define $ L_{\mu}^A = ( \partial_{A_\mu} H)|_{A=0} $ and so on. 


> However, shouldn't this second line just be 
$$ H= H_0 + L_\mu^A A_\mu + 1/2 L_{\mu v}^{AA} A_\mu A_v +... $$

Both expressions look wrong to me. Theirs for the reason you already know, yours because you still have to integrate over the dummy position $\vec r'$ in $A(\vec r',t)$. So unless there is some implicit convolution in your notation, yours looks wrong too.

The expansion is obtained by using the Taylor series expansion of the exponential:
$$
H = \sum_{i,j} \sum_{a,b} t_{a,b}   \left(1 
+ i \int_{\textbf{R}_{j,b}}^{\textbf{R}_{i,a}} dr'_{\mu} A_{\mu} (\textbf{r}',t ) 
-
\frac{1}{2}\int_{\textbf{R}_{j,b}}^{\textbf{R}_{i,a}} dr'_{\mu} A_{\mu} (\textbf{r}',t )
\int_{\textbf{R}_{j,b}}^{\textbf{R}_{i,a}} dr''_{\mu} A_{\mu} (\textbf{r}'',t ) +\ldots
\right) c_{i,a}^\dagger c_{j,b}
$$