Sorry I didn’t see this earlier. My memory is vague, and probably colored by subsequent events and results, but here’s how I recall things happening.
Since I had read and enjoyed Lazard’s paper on one-dimensional formal group (laws), which dealt with the case of a base field of characteristic $p$, I decided to look at formal groups over $p$-adic rings. For whatever reason, I wanted to know about the endomorphism rings of these things, and gradually recognized the similarity between, on the one hand, the case of elliptic curves and their supersingular reduction mod $p$, when that phenomenon did occur, and, on the other hand, formal groups over $p$-adic integer rings of higher height than $1$.
I had taken, or sat in on, Tate’s first course on Arithmetic on Elliptic Curves, and was primed for all of this. In addition, I was aware of Weierstrass Preparation, and the power it gave to anyone who wielded it. And in the attempt to prove a certain result for my thesis, I had thought of looking at the torsion points on a formal group, and I suppose it was clear to me that they formed a module over the endomorphism ring. Please note that it was not my idea at all to use them as a representation module for the Galois group.
But Tate was looking over my shoulder at all times, and no doubt he saw all sorts of things that I was not considering. At the time of submission of my thesis, I did not have a construction of formal groups of height $h$ with endomorphism ring $\mathfrak o$ equal to the integers of a local field $k$ of degree $h$ over $\Bbb Q_p$. Only for the unramified case, and I used extremely tiresome degree-by-degree methods based on the techniques of Lazard. Some while after my thesis, I was on a bus from Brunswick to Boston, and found not only that I could construct formal groups in all cases that had this maximal endomorphism structure, but that one of them could take the polynomial form $\pi x+x^q$. Tate told me that when he saw this, Everything Fell Into Place. The result was the wonderful and beautiful first Lemma in our paper, for which I can claim absolutely no responsibility. My recollection, always undependable, is that the rest of the paper came together fairly rapidly. Remember that Tate was already a master of all aspects of Class Field Theory. But if the endomorphism ring of your formal group is $\mathfrak o$ and the Tate module of the formal group is a rank-one module over this endomorphism ring, can the isomorphism between the Galois group of $k(F[p^\infty]])$ over $k$ and the subgroup $\mathfrak o^*\subset k^*$ fail to make you think of the reciprocity map?
Given $\phi(x)\in\mathbb{R}[[x]]$, with $\phi(0)=1$, we have defined $g(x):=\int^x_0{dt\over \phi(t)}$, $f:=g^{-1}$ and $$F(x,y)=f\big(g(x)+g(y)\big)=\sum_{n=0}^\infty \psi_n(x) {y^n\over n!}\in\mathbb{R}[[x,y]].$$ Let's write a recursion for the coefficient sequence $\psi_n=\partial_y^nF(x,0)\in\mathbb{R}[[x]]$, solving by series the differential equation satisfied by $F$,
$$\cases{\phi(x)\, F_x(x,y)=\phi(F(x,y))\\ F(x,0)=x\ .}$$
One finds $\psi_0=x$, $\psi_1=\phi,\dots$ . Let's take $\partial_y^n$ at $ {y=0}$ on both sides. Faà di Bruno:
$$\partial_y^n\big( \phi\circ F\big)\big|_{y=0}=\Big(\sum_{\alpha\in\operatorname{par}[n]} \phi_y^{(|\alpha|)} (F)\, \prod_{s\in\alpha} \partial_y^{|s|}F \Big) \ \Big|_{y=0} =\sum_{\alpha\in\operatorname{par}[n]} \phi^{(|\alpha|)}(x) \prod_{s\in\alpha} \psi_{|s|}(x) ,$$
(Legenda: The sum is indexed on the set of all partitions of $[n]:=\{1,2,\dots,n\}$, and $|\cdot|$ denotes cardinality. The latter equality comes from $F(x,0)=x$ and $\partial_y^{j}F(x,y)\big|_{y=0}=\psi_j$).
Now we isolate the term $\phi'\psi_n$, that corresponds to the partition $\alpha$ into a single class, from the terms of the sum indexed on the set of non-trivial partitions, with $|\alpha|>1$, denoted $\operatorname{par}^*[n]$. Note that each of these terms contains more than one factor $\psi_j$.
$$\phi\psi'_n -\phi' \psi_n =\sum_{\alpha\in\operatorname{par}^*[n]} \phi^{(|\alpha|)} \prod_{s\in\alpha} \psi_{|s|} .$$
Multiplying by the integrating factor $\phi^{-2}$ , and since $\psi_n(0)=0$, for $n>1$
$$ \psi_n(x) =\phi(x)\int_0^x\big(\!\sum_{\alpha\in\operatorname{par}^*[n]} \phi^{(|\alpha|)} \prod_{s\in\alpha} \psi_{|s|}\,\big)\phi^{-2}\ dt .$$
It is now clear by complete induction that for any $n\ge1$, $\psi_n$ is equal to $\phi$ times a series with positive coefficients, proving your conjecture.
Best Answer
Okay, here's a few words about the relation between the $L$-series and the formal group. In general, if $F(X,Y)$ is the formal group law for $\hat G$, then there is an associated formal invariant differential $\omega(T)=P(T)dT$ given by $P(T)=F_X(0,T)^{-1}$. Formally integrating the power series $\omega(T)$ gives the formal logarithm $\ell(T)=\int_0^T\omega(T)$. The logarithm maps $\hat G$ to the additive formal group, so we can recover the formal group as $$ F(X,Y) = \ell^{-1}(\ell(X)+\ell(Y))$$. (See, e.g., Chapter IV of Arithmetic of Elliptic Curves for details.)
Now let $E$ be an elliptic curve and $\omega=dx/(2y+a_1x+a_3)$ be an invariant differential on $E$. If $E$ is modular, say corresponding to the cusp form $g(q)$, then we have (maybe up to a constant scaling factor) $\omega = g(q) dq/q = \sum_{n=1}^{\infty} a_nq^{n-1}$. Eichler-Shimura tell us that the coefficients of $g(q)$ are the coefficients of the $L$-series $L(s)=\sum_{n=1}^\infty a_n n^{-s}$. Integrating $\omega$ gives the elliptic logarithm, which is the function you denoted by $f$, i.e., $f(q)=\sum_{n=1}^\infty a_nq^n/n$, and then the formal group law on $E$ is $F(X,Y)=f^{-1}(f(X)+f(Y))$.
To me, the amazing thing here is that the Mellin transform of the invariant differential gives the $L$-series. Going from the invariant differential to the formal group law via the logarithm is quite natural.