I suppose it's bad form to answer one's own MO question, but I now have an almost complete solution to this one. I can prove:
1.----H is always in the ring S generated by the [j].
2.----The same holds for G except perhaps when l=15 mod 16. (In "More questions involving characteristic 2 theta series identities" I provide some experimental evidence when l=15 mod 16.)
To prove 1. note that I gave a formula in my question expressing H as a quotient of elements of S. Now I have made a study of the variety V consisting of the zeros of the polynomial relations between the various [j]. V is a curve; when l>3 it has exactly l+1 singular points, each of which is an ordinary multiple point of multiplicity (l-1)/2. Using my formula for H, I can show that it has ord at least 0 at every non-singular point of V, and ord> 0 at every branch centered at every singular point. So it lies in all the local rings of S.
EDIT:NOT SO--the condition of being in the local ring at a singular point is more stringent.
For a correct argument see the FINAL EDIT below.
To prove 2. let C be the sum of the x^(ln) where n runs over all (non-zero) integers of the form (square) or 2(square) or l(square) or 2l(square). Note that C^2+C is G+H. So in view of 1. it suffices to show that C is in S. In my previous answer I indicated why this is true when l=1 mod 4 or l=3 mod 8, writing C explicitly as a polynomial in the [j]. I will edit this answer shortly to handle the more difficult case l=7 mod 16.
EDIT: Suppose now l=7 mod 16. Here's a proof that C lies in S. Let T, contained in a product of 4 copies of Z/l, consist of all (r1,r2,r3,r4) other than (0,0,0,0) with (r1^2)+(r2^2)+(r3^2)+(r4^2)=0. There is a group of order (24)(16)=384 acting on T by permutation of co-ordinates and sign changes of co-ordinates. Using the fact that l=7 mod 8, we find that every orbit has size 384 or 192 or 64. Call an orbit "small" if it has size 192 or 64. We shall show that C is a sum of terms attached to the small orbits. To each small orbit attach the power series [r1][r2][r3][r4] where (r1,r2,r3,r4) is an orbit representative; this is independent of the representative. I'll show that C is the sum of these contributions. Clearly every exponent appearing in the sum of these power series is
divisible by l. It remains to show that x^ln appears in the sum if and only if n is the product of a non-zero square by 1,2,l or 2l.
Now the coefficient of x^ln in [r1][r2][r3][r4] is the mod 2 reduction of M where M is the number of integer 4-tuples (a,b,c,d) satisfying:
(1)---(a^2)+(b^2)+(c^2)+(d^2)=0
(2)---(a,b,c,d) reduces to (r1,r2,r3,r4) mod l
Modulo 2, M is (1/64)(the number of (a,b,c,d) satisfying (1) and reducing to an element in the orbit of (r1,r2,r3,r4)). So the sum of the M, modulo 2, is the number of (a,b,c,d) satisfying (1) and reducing to a point in some small orbit. Also the number of (a,b,c,d)
satisfying (1) and reducing to a point in an orbit of size 384 is clearly a multiple of 384.
So the coefficient of x^ln in our sum is the mod 2 reduction of (1/64)(the number of (a,b,c,d) that satisfy (1) and do not reduce to (0,0,0,0)).
Let R(n) be the number of representations of n as a sum of 4 squares. We have just shown that the coefficient of x^ln in our sum is the mod 2 reduction of (1/64)(R(ln)-R(n/l)). It remains to show that (1/64)(R(ln)-R(n/l)) is odd precisely when n is the product of a square by 1,2,l or 2l. Jacobi proved that R(n) is 8(the sum of the divisors of n) when n is odd, and 24(the sum of the odd divisors) when n is even. So (1/8)(R(ln)-R(n/l)) is a product of
local factors, one from each prime. The factor attached to 2 is 1 or 3. That attached to l is l^t(1+l) where t=(ord_l)(n). Since l=7 mod 16, this is 8(an odd number). Finally that attached to an odd prime p other than l is the sum of the divisors of p^s where s is (ord_p)(n). This factor is odd just when s is even, and the result follows.
A couple of remarks. When l=15 mod 16 the same argument shows that the sum we've constructed is not C, but 0. Also an orbit is small precisely when it has a representative with r1=r2 or a representative with r4=0.
FINAL EDIT: I now have an answer I'm prepared to accept, unless some spoilsport finds a flaw; it shows that G,H (and F) all lie in the subring S of Z/2[[x]] generated by the [j]
irrespective of l. Unlike the approach taken in the last edit which exhibited G+H explicitly as a polynomial in the [j], (except when l is 15 mod 16), this one doesn't seem to give nice explicit formulas. I'll be using results from other MO questions of mine, and some further results in manuscript. Let K be an algebraic closure of Z/2, and S' be the subring of K[[x]] generated over K by the [j]. It,s enough to show that G,H and F lie in S'.
First I show that they're all in the field of fractions, L, of S'. In another MO post I wrote H as a quotient of 2 elements of S. To handle F I use the following:
(1)___For l>3, Spec(S') is a curve with l+1 singular points, among them the maximal ideal m generated by [1],...,[l-1]. These are ordinary singular points of multiplicity (l-1)/2.
(2)___There is a group of automorphisms of S'/K isomorphic to PSL_2(Z/l). These automorphisms stabilize the space spanned by [0],...,[l-1] and act transitively on the (l+1)(l-1)/2 valuation rings in L/K containing the local rings at the singular points. The group is generated by the maps [j]-->[rj], r prime to l, [j]-->a^(j^2) [j] where a is an l'th
root of unity in L, and a sort of characteristic 2 "Fourier transform".
Now the maps [j]-->[rj] and [j]-->a^(j^2) [j] generate a subgroup B of PSL_2 of order l(l-1)/2, and my "quotient formula for H" shows that B fixes H. So the orbit of H under PSL_2 has size at most l+1. A rather formal calculation with the "Fourier transform" shows that the orbit consists of H and the F(ax) where a^l=1. I claim that each of these elements lies in the local ring of m on S'. For H this is easy; H has ord l^2 at each valuation ring containing m. Taking E to be the sum of [1],...[(l-1)/2] we find that E+E^4=F+H. So F is in this local ring as well, and the result follows easily for each F(ax). The fact that PSL_2
acts transitively on the singular points now shows that H and the F(ax) lie in the local ring at every singular point. Also the quotient formula for H shows that H has ord 0 at every non-singular point, and the same then holds for the F(ax). Thus H and the F(ax) are in S'; this corrects the argument I gave earlier.
I now turn to G. There is a degree l+1 2-variable symmetric polynomial P over Z/2 with P(F,G)=0. Furthermore P(z,G) is monic of degree l+1, and has H and the F(ax) as roots. Also the constant term of P(z,G) is G^(l+1), while the coefficient of z is G+ higher degree terms. Since the product of H and the F(ax), as well as the l'th symmetric function of H and the F(ax), are in S', both G^(l-1) and G+... are in S'. Now over K these 2 elements generate a field between K(G^(l+1)) and K(G); since G+... is in this field it is all of K(G), and G is in L. Also G^(l+1), as the product of H and the F(ax), is fixed by PSL_2.
Since every homomorphism from PSL_2 to the l+1 th roots of unity is trivial, G is fixed by
PSL_2.
At the valuation rings lying over m, G has ord l. So G is in the local ring of m, and consequently in the local ring at every singular point. Furthermore, like H and the F(ax), G
has ord 0 at the non-singular points. So it is in S'. (Note also that like H and the F(ax), G has poles of order 12 at every valuation ring in L/K that doesn't contain S').