The functionality of `EtaProduct`

seems to be restricted to creating modular functions.
For the purposes of the posted question there is no implemented support.
(At least as far as my quick search failed to find any.)

In the following, there will be some lines of code, written and dedicated for the case
$$
S_6(\ \Gamma_0(24)\ )\ .
$$
(The question was a good occasion for me to learn about the subject,
i had to do something in a neighbour room, and was even extermely curious to see the results
for $S_6(\Gamma_0(24))$.)

(The code is written,
so that it could be reused, at least theoretically.)

As a background for the implementation, it is simplest for me to use the setting (definitions, notations, and results)
from
[Günther Köhler, Eta Products and Theta Series Identities, Springer Monographs in Mathematics, Springer, 2011],
in the sequel
[GK2011], for short.
For the convenience of the reader, i will try to also illustrate the intermediate steps by the structure found so far.

ONE.

Chapter 2 in [GK2011], conditions for holomorphic eta products, is telling us which data
can be used to build relevant eta products.

It is convenient to consider an eta product associated to a symbolic sum of the shape
$$
x = \sum_{d|24}x(d)[d]\in{\mathbb Z}[D]
$$
in the free abelian group, i.e. free ${\mathbb Z}$--module
with the basis $[d]$, $d$ being a pure symbol running in the set $D$ of divisors of $24$.
(There will be no $N=24$, for me it is better to fix ideas and use the $24$ explicitly.)

(The additive notation is better suited, the multiplicative notation
$\prod_{d|24}[d]^{x(d)}$
represents small the important data $(x(d))_d$, as an upper index.)

Associated to $x$ there is the following eta product:
$$
\eta(x)(z)=\prod{d|24}\eta(\; dz\;)^{x(d)}\ .
$$
The conditions of holomorphicity at the cusps (of the form $1/c$, where $c$ divides $24$,) can be restated as
$$
24\sum_{d|24}\frac{\text{gcd}(c,d)^2}{cd}\cdot x(d)\ge 0\ ,
$$
for all $c$ with $c|24$.
(It is
This isolates a symmetric matrix with entries
$$
24\frac{\text{gcd}(c,d)^2}{cd}\ ,\ c,d\in D\ ,
$$
which are integers.

We fix an order on $D$. In the following code, it corresponds to
the entries of the $8\times 1$ matrix `d24`

. (The number $24$ has `euler_phi(24)`

$=8$ divisors.)

The above entries are then arranged corresponding to this order and build a matrix $B=$`B24`

.

We search column vectors $x$ (with entries $x(d)$) so that the holomorphicity conditions are satisfied:
$$
Bx\ge 0\ .
$$
(The $0$ is the $8\times 1$ zero matrix and the inequality is considered componentwise.)
This isolates a cone in ${\mathbb Z}^8$.

The eta function behaves like a modular form of weight $1/2$ for some multiplier system
that can be technically described. So $\eta(x)$ behaves like a modular form of weight
$$
\frac 12 \sum_{d|24} x(d)\ .
$$
In our case, we are searching modular (and even cusp) forms of weight $6$, so the condition is
$$
12 = \sum_{d|24} x(d)\ .
$$
We intersect the above cone $Bx\ge 0$ with this hyperplane.
The resulted simplex has "known extremal points" $v$, that can be explicitly described
for "the same problem", when we replace $24$ by its highest prime power divisors
$8$, and respectively $3$. Such extremal points $v$ build correspondingly (quadratic) matrices `V8`

and `V3`

with `euler_phi(8)`

$=4$, and respectively `euler_phi(3)`

$=2$ rows and columns.
The following routine `dBV`

computes the data $d,B,V$ (order of the divisors,
matrix $B$ describing the cone simplex $Bx\ge 0$, and some scaling to integers of the extremal points of the cone $V$).
We will soon have:

```
sage: dBV(8)
(
[1] [8 4 2 1] [ 2 -2 0 0]
[2] [4 8 4 2] [-1 5 -2 0]
[4] [2 4 8 4] [ 0 -2 5 -1]
[8], [1 2 4 8], [ 0 0 -2 2]
)
sage: dBV(3)
(
[1] [3 1] [ 3/2 -1/2]
[3], [1 3], [-1/2 3/2]
)
```

The third matrix is in both cases so that the sum of entries in each column is one.

Then the data for $24$ is collected by building the tensor product of the above matrices.
The results are block matrices:

```
sage: dBV(24)
(
[ 1] [24 8|12 4| 6 2| 3 1] [ 3 -1| -3 1| 0 0| 0 0]
[ 3] [ 8 24| 4 12| 2 6| 1 3] [ -1 3| 1 -3| 0 0| 0 0]
[--] [-----+-----+-----+-----] [---------+---------+---------+---------]
[ 2] [12 4|24 8|12 4| 6 2] [-3/2 1/2|15/2 -5/2| -3 1| 0 0]
[ 6] [ 4 12| 8 24| 4 12| 2 6] [ 1/2 -3/2|-5/2 15/2| 1 -3| 0 0]
[--] [-----+-----+-----+-----] [---------+---------+---------+---------]
[ 4] [ 6 2|12 4|24 8|12 4] [ 0 0| -3 1|15/2 -5/2|-3/2 1/2]
[12] [ 2 6| 4 12| 8 24| 4 12] [ 0 0| 1 -3|-5/2 15/2| 1/2 -3/2]
[--] [-----+-----+-----+-----] [---------+---------+---------+---------]
[ 8] [ 3 1| 6 2|12 4|24 8] [ 0 0| 0 0| -3 1| 3 -1]
[24], [ 1 3| 2 6| 4 12| 8 24], [ 0 0| 0 0| 1 -3| -1 3]
)
```

(The print shows `d24`

, the chosen order of the divisors of $24$, the matrix
`B24`

, and the columns of the next block matrix `V24`

are extremal points for the cone.)
The product of the last two matrices is:

```
sage: d24, B24, V24 = dBV(24)
sage: B24 * V24
[48 0 0 0 0 0 0 0]
[ 0 48 0 0 0 0 0 0]
[ 0 0 96 0 0 0 0 0]
[ 0 0 0 96 0 0 0 0]
[ 0 0 0 0 96 0 0 0]
[ 0 0 0 0 0 96 0 0]
[ 0 0 0 0 0 0 48 0]
[ 0 0 0 0 0 0 0 48]
```

An eta product, leading to a (holomorphic) modular form of weight $6$
is thus obtained by taking a linear combination with coefficients (building a vector) $t\ge 0$
of the columns of `V24`

, such that the scalars used in it sum to $12=2\cdot 6$.
(Since each column of `V24`

has entries adding to one.)

We will get an $x$ of the shape $x=Vt$.

Let us type some code first, that collects the matrices $d, B,V$, this will fix (for me) the frame:

```
def dBV( N ):
"""deliver three matrices, d, B, and V,
d records the divisors of N, computed in a specific order,
that will reflect the same order used for B, V.
B implements the / a positivity condition "pos",
a solution of "pos" is a column x with Bx >= 0,
and V is the list of extremal points (as columns) for "pos" whith rational coordinates,
and such that each column sum is one.
This will help us to construct eta products as follows.
The components of d, considered as eta argument multiplicators,
and the corresponding components of a solution x
(x has components in ZZ), considered as corresponding eta powers,
describe a relevant holomorphic (at cusps) eta product.
It is a cusp form, iff Bx > 0.
"""
if N.is_prime_power():
p, ppower = N.is_prime_power( get_data=True ) # so N = p^ppower
nv = ppower+1
d = matrix( QQ, nv, 1 )
b = matrix( QQ, nv, nv )
v = matrix( QQ, nv, nv )
for j in range(nv):
d[j, 0] = p^j
for k in range(nv):
b[j, k] = p^( ppower + 2*min(j,k) - j - k )
for j in range(nv):
if j == 0:
v[0, 0] = p*(p-1)
v[1, 0] = 1-p
elif j == nv-1:
v[j-1, j] = 1-p
v[j , j] = p*(p-1)
else:
v[j-1, j] = -p
v[j , j] = 1+p^2
v[j+1, j] = -p
return d, b, (1/(p-1)^2) * v
else:
prime_powers = [ p^r for p, r in N.factor() ]
d = matrix( QQ, 1, 1, [1,] )
v = matrix( QQ, 1, 1, [1,] )
b = matrix( QQ, 1, 1, [1,] )
for q in prime_powers:
dq, bq, vq = dBV( q )
d = d.tensor_product( dq )
b = b.tensor_product( bq )
v = v.tensor_product( vq )
return d, b, v
d24, B24, V24 = dBV( 24 )
dd24 = matrix( QQ, 8, 1, [ 24/d24[k, 0] for k in range(8) ] )
print "d24 is:\n%s\n" % d24
print "dd24 is:\n%s\n" % dd24
print "B24 is:\n%s\n" % B24
print "V24 is:\n%s\n" % V24
print "B24 * V24 is:\n%s\n" % ( B24*V24 )
```

The above calculates, as already mentioned:

```
d24 is:
[ 1]
[ 3]
[--]
[ 2]
[ 6]
[--]
[ 4]
[12]
[--]
[ 8]
[24]
dd24 is:
[24]
[ 8]
[12]
[ 4]
[ 6]
[ 2]
[ 3]
[ 1]
B24 is:
[24 8|12 4| 6 2| 3 1]
[ 8 24| 4 12| 2 6| 1 3]
[-----+-----+-----+-----]
[12 4|24 8|12 4| 6 2]
[ 4 12| 8 24| 4 12| 2 6]
[-----+-----+-----+-----]
[ 6 2|12 4|24 8|12 4]
[ 2 6| 4 12| 8 24| 4 12]
[-----+-----+-----+-----]
[ 3 1| 6 2|12 4|24 8]
[ 1 3| 2 6| 4 12| 8 24]
V24 is:
[ 3 -1| -3 1| 0 0| 0 0]
[ -1 3| 1 -3| 0 0| 0 0]
[---------+---------+---------+---------]
[-3/2 1/2|15/2 -5/2| -3 1| 0 0]
[ 1/2 -3/2|-5/2 15/2| 1 -3| 0 0]
[---------+---------+---------+---------]
[ 0 0| -3 1|15/2 -5/2|-3/2 1/2]
[ 0 0| 1 -3|-5/2 15/2| 1/2 -3/2]
[---------+---------+---------+---------]
[ 0 0| 0 0| -3 1| 3 -1]
[ 0 0| 0 0| 1 -3| -1 3]
B24 * V24 is:
[48 0 0 0 0 0 0 0]
[ 0 48 0 0 0 0 0 0]
[ 0 0 96 0 0 0 0 0]
[ 0 0 0 96 0 0 0 0]
[ 0 0 0 0 96 0 0 0]
[ 0 0 0 0 0 96 0 0]
[ 0 0 0 0 0 0 48 0]
[ 0 0 0 0 0 0 0 48]
```

TWO:

We want moreover a cusp form resulting from the construction, so for $Bx>0$, $x=Vt$,
we need correspondingly $t>0$.
In order to get some runnable code, some examples, i will start with the Ansatz,
that t has natural entries $\ge 1$. (So there will be no component $1/2$ or $1/3$ or...)
This facilitates the enumeration. Then the eight entries of $t$ are at least one each, they sum to $12$,
so we have four further points to distribute among the components.

The following simple code collects the candidates:

```
etaDataList = []
t0 = matrix( QQ, 8, 1, [ 1 for _ in range(8) ] )
tList = []
C = [ [c0, c1, c2, c3]
for c0 in [ 0..7 ]
for c1 in [ c0..7 ]
for c2 in [ c1..7 ]
for c3 in [ c2..7 ] ]
for c in C:
t = t0 + matrix( QQ, 8, 1, [ c.count( _ ) for _ in range(8) ] )
x = V24 * t
if not all( [ entry in ZZ for entry in ( d24.transpose() * t ) . list() ] ): continue
# if not all( [ entry in ZZ for entry in ( dd24.transpose() * t ) . list() ] ): continue
if not all( [ entry in ZZ for entry in ( x ) . list() ] ): continue
tList.append( t )
etaDataList.append( ( [ ZZ(entry) for entry in x.list() ], d24 ) )
```

(Decommenting the commented condition leads to the same list.)
Here, `len( etaDataList )`

and `len( tList )`

are both $170$.
So we have $170$ candidates to play with.

(We do not need the `tList`

, and also need from `etaDataList`

only the data comming from the $x$-component.
But for $170$ entries this is not expensive.)

THREE:

We are in the position to construct some eta products.
Let us do this with one entry in the collected `etaDataList`

, for instance with

```
ed = etaDataList[ 11 ]
```

This is

```
sage: ed
(
[ 1]
[ 3]
[--]
[ 2]
[ 6]
[--]
[ 4]
[12]
[--]
[ 8]
[6, -2, -3, 1, 7, 7, -2, -2], [24]
)
```

Associated to it we have the linear combination
$$ \sum x(d)[d] =
6[1] - 2[3] - 3[2] + [6] + 7[4] + 7[12] - 2[8] - 2[24]
\in{\mathbb Z}[D]\ .
$$
The eta product for it is
$$ \eta(x)(z) =
\eta(\; z\; )^6
\eta(\; 3z\; )^{-2}
\eta(\; 2z\; )^{-3}
\eta(\; 6z\; )
\eta(\; 4z\; )^7
\eta(\; 12z\; )^7
\eta(\; 8z\; )^{-2}
\eta(\; 24z\; )^{-2}\ .
$$
The corresponding $q$-product representation, $q = \exp( \; 2\pi i\; z\; )$,
has $d\cdot z\leftrightarrow q^d$, and $\eta(z)\leftrightarrow q^{1/24}f(q)$,
where $f$ "is" `f`

from

```
PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )
```

which gives

```
sage: f
1 - q - q^2 + q^5 + q^7 - q^12 - q^15 + q^22 + q^26 - q^35 - q^40 + q^51 + q^57 - q^70 - q^77 + q^92 + O(q^100)
```

and $\eta(d\cdot z)\leftrightarrow q^{d/24}f(q^d)$,
and in the eta product we have to cumulate the power-of-$q$-prefactor.
To the above element $x=6[1] - 2[3] - 3[2] + [6] + 7[4] + 7[12] - 2[8] - 2[24]
\in{\mathbb Z}[D]$ there corresponds...

```
ed = etaDataList[ 11 ]
x, d = ed
PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )
fx = PSR( prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) \
* \
q^ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 ) )
fx
```

This gives a "long" $q$-expansion, and we can / have to truncate it here to fit somehow:

```
sage: fx + O(q^10)
q^2 - 6*q^3 + 12*q^4 - 6*q^5 - 13*q^6 + 42*q^7 - 72*q^8 + 30*q^9 + O(q^10)
```

"Is" this element in the space of cusp forms $S_6(\Gamma_0(24))$, at least
if we test this for the first $100$ coefficients in the $q$-expansion?!
For this, we associate the space `CF = CuspForms( Gamma0(24), 6 )`

of cusp forms,
each cusp form (in the basis of this space) is mapped to its $q$--expansion,
truncated modulo $O(q^100)$, then mapped to $\mathbb Q^{100}$ by picking the coefficients
with respect to $1,q,q^2,\dots, q^{99}$ modulo $O(q^{100})$ .

This is a vector space $W=$`WCF`

of dimension $16$, let us construct it:

```
CF = CuspForms( Gamma0(24), 6 )
# CF is:
# Cuspidal subspace of dimension 16
# of Modular Forms space of dimension 24
# for Congruence Subgroup Gamma0(24)
# of weight 6 over Rational Field
def fit_to_length( myList, length ):
lenMyList = len( myList )
if lenMyList == length:
return myList
elif lenMyList < length:
return myList + ( length - lenMyList ) * [0,]
else:
return myList[:length]
VS = VectorSpace( QQ, 100 )
W = VS.subspace_with_basis( [ VS( fit_to_length( cf.qexp(100).list(), 100 ) )
for cf in CF.basis() ] )
```

Then the coordinates of the eta product $\eta(x)$ in the sage basis of `CF`

are as follows:

```
Wf = VS( fit_to_length( fx.list(), 100 ) )
coord = W.coordinate_vector( Wf )
print coord
```

This gives:

```
(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)
```

In particular, the constructed eta product lives in the space of cusp forms `CF`

$=S_6(\; \Gamma_0(24)\; )$.

FOUR:

Let us print the first coefficients, say working modulo $O(q^{21})$ in the `CF`

basis and in the eta
product $\eta(x)$. The $16$ elements in the basis are (starting with)

```
for cf in CF.basis():
print cf.qexp( 21 )
```

This gives:

```
q + 48*q^19 + O(q^21)
q^2 + 39*q^18 + O(q^21)
q^3 + 30*q^19 + O(q^21)
q^4 + 4*q^16 - 30*q^20 + O(q^21)
q^5 - 20*q^19 + O(q^21)
q^6 - 12*q^18 + O(q^21)
q^7 - 7*q^19 + O(q^21)
q^8 - 6*q^16 + O(q^21)
q^9 - 20*q^19 + O(q^21)
q^10 - 15*q^18 + O(q^21)
q^11 - 11*q^19 + O(q^21)
q^12 - 4*q^16 + 4*q^20 + O(q^21)
q^13 - 7*q^19 + O(q^21)
q^14 - 6*q^18 + O(q^21)
q^15 - 5*q^19 + O(q^21)
q^17 - q^19 + O(q^21)
```

And of course, from the $q$--expansion of $\eta(x)$,

```
sage: fx + O(q^21)
q^2 - 6*q^3 + 12*q^4 - 6*q^5 - 13*q^6 + 42*q^7 - 72*q^8 + 30*q^9 + 70*q^10
- 102*q^11 + 108*q^12 - 12*q^13 - 200*q^14 + 114*q^15 + 48*q^16 - 108*q^17
+ 345*q^18 - 210*q^19 + 72*q^20 + O(q^21)
```

(Output was manually rearranged.)

The coefficients in

```
(0, 1, -6, 12, -6, -13, 42, -72, 30, 70, -102, 108, -12, -200, 114, -108)
```

are exactly the $16$ coefficients of the $q$--expansion in degrees
$1,2,\dots, 14,15,\ 17$. (Which are the "(strating) orders" for the basis elements
of `CF`

.)

FIVE:

Let us see that all $170$ eta products found in TWO (variable
`etaDataList`

)
are also in `CF`

, at least regarding the test of the first
$100$ coefficients.

```
PSR.<q> = PowerSeriesRing( QQ, default_prec=100 )
f = qexp_eta( PSR, 100 )
VS = VectorSpace( QQ, 100 )
W = VS.subspace_with_basis( [ VS( fit_to_length( cf.qexp(100).list(), 100 ) )
for cf in CF.basis() ] )
count = 0
for ed in etaDataList:
count += 1
x, d = ed
qPower = ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 )
fx = PSR( q^qPower * prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) )
Wfx = VS( fit_to_length( fx.list(), 100 ) )
coord = W.coordinate_vector( Wfx )
print "%3s x=%-30s -> %s" % ( count, str(x), str(coord) )
```

This gives a long list of coordinates for the found eta products
with respect to the standard base of `CF`

...

The first ten entries, reshaped in the print are:

```
count = 0
for ed in etaDataList[:10]:
count += 1
x, d = ed
qPower = ZZ( sum( [ d[k, 0]*x[k] for k in range(8) ] ) / 24 )
fx = PSR( q^qPower * prod( [ f( q^ZZ(d[k, 0]) )^x[k] for k in range(8) ] ) )
Wfx = VS( fit_to_length( fx.list(), 100 ) )
coord = W.coordinate_vector( Wfx )
print "%3s x=%-30s\n %s" % ( count, str(x), str(coord) )
```

This gives:

```
1 x=[12, -4, -4, 4, 2, 2, 0, 0]
(0, 1, -12, 58, -132, 81, 276, -668, 468, 294, -876, 810, 24, -904, 324, -936)
2 x=[8, 0, -2, 2, 2, 2, 0, 0]
(0, 1, -8, 22, -16, -27, 40, -20, 96, -138, -8, 54, -128, 392, -216, -192)
3 x=[6, -2, 5, 1, -1, 3, 0, 0]
(0, 1, -6, 4, 42, -81, -78, 304, -90, -354, 426, -324, -204, 1040, -486, 180)
4 x=[10, -6, -5, 11, 3, -1, 0, 0]
(0, 1, -10, 40, -74, 27, 158, -344, 282, 78, -442, 432, -52, -256, 54, -564)
5 x=[4, 4, 0, 0, 2, 2, 0, 0]
(0, 1, -4, 2, 4, 9, -4, -44, 12, 6, 92, 18, -88, -40, -180, -24)
6 x=[2, 2, 7, -1, -1, 3, 0, 0]
(0, 1, -2, -8, 14, 27, -26, -56, -30, 78, 142, 0, -68, -256, -162, 60)
7 x=[6, -2, -3, 9, 3, -1, 0, 0]
(0, 1, -6, 12, -6, -9, 18, -32, 54, -66, 42, 36, -108, 176, -198, -108)
8 x=[0, 0, 14, -2, -4, 4, 0, 0]
(0, 1, 0, -14, 0, 81, 0, -236, 0, 294, 0, 162, 0, -904, 0, 0)
9 x=[4, -4, 4, 8, 0, 0, 0, 0]
(0, 1, -4, -2, 28, -27, -52, 124, -60, -138, 284, -162, -136, 392, -324, 120)
10 x=[8, -8, -6, 18, 4, -4, 0, 0]
(0, 1, -8, 26, -40, 9, 88, -188, 168, 6, -200, 234, -80, -40, -72, -336)
```