This notebook contains additional material to the Mathematica GuideBooks  (Springer-Verlag, New York, 2004) http://www.MathematicaGuideBooks.org . It uses the stylesheet of the Mathematica GuideBooks, which is embedded.

Symbolics Volume, Chapter 3

The elliptic theta functions (Subscript[ϑ, a])(z,q) obey many identities. For explicit values of z, typical examples are:

(Subscript[ϑ, 3])((π/4),q)-(Subscript[ϑ, 4])(0,(q^4))=0

2 (Subscript[ϑ, 1])'(0,(q^4))+(Subscript[ϑ, 3])'((π/4),q)=0

(Subscript[ϑ, 1])(-(π/2),q) (Subscript[ϑ, 1])(-(π/4),q)+2 (Subscript[ϑ, 1])(-(π/4),(q^2)) (Subscript[ϑ, 3])((π/4),(q^2))=0

2 (Subscript[ϑ, 3])(-(π/2),(q^2)) (Subscript[ϑ, 3])'(-(π/4),q)-(Sqrt[2]) (Subscript[ϑ, 1])(-(π/2),q) (((Subscript[ϑ, 1])(-(π/4),q))^3)=0

(Sqrt[2]) (Subscript[ϑ, 1])(-(π/2),q) (Subscript[ϑ, 1])(-(π/4),q) (Subscript[ϑ, 3])(-(π/4),q)-2 (Subscript[ϑ, 1])'(0,(q^2))=0

(Sqrt[2]) ((((Subscript[ϑ, 1])(-(π/2),q))^2)+2 (((Subscript[ϑ, 1])(-(π/4),q))^2)) (Subscript[ϑ, 1])(-(π/4),(q^2))-2 (Subscript[ϑ, 1])(-(π/2),(q^4)) ((((Subscript[ϑ, 3])((π/8),q)+(Subscript[ϑ, 4])((π/8),q)))^2)=0

(Subscript[ϑ, 1])'(0,q) (((Subscript[ϑ, 1])((π/3),(q^2)))^2)+(Subscript[ϑ, 1])(-(π/2),q) (((Subscript[ϑ, 1])((π/3),q))^2) (((Subscript[ϑ, 1])((π/6),q))^2)=0

(((Subscript[ϑ, 1])(-(π/3),(q^2)))^2) (((Subscript[ϑ, 1])(-(π/3),q))^4)-9 (((Subscript[ϑ, 1])'(-(π/6),(q^4)))^2)+9 (((Subscript[ϑ, 3])'(-(π/3),(q^4)))^2)=0

(((Subscript[ϑ, 1])(-(π/3),q))^2) (((Subscript[ϑ, 1])(-(π/2),q))^6)+8 (((Subscript[ϑ, 1])(-(π/3),q))^2) (((Subscript[ϑ, 1])(-(π/6),q))^6)+8 (((Subscript[ϑ, 1])(-(π/3),q))^6) (((Subscript[ϑ, 1])(-(π/6),q))^2)-144 (((Subscript[ϑ, 1])(-(π/3),(q^2)))^2) (((Subscript[ϑ, 1]^′)(-(π/3),(q^2)))^2)=0.

20 (((Subscript[ϑ, 1])(-(π/3),q))^4) (((Subscript[ϑ, 1])'(0,(q^2)))^2)-3 (Subscript[ϑ, 1])'(-(π/3),q) (((Subscript[ϑ, 1])(-(π/2),q))^7)-72 (((Subscript[ϑ, 1])'(-(π/3),(q^2)))^2) (((Subscript[ϑ, 1])(-(π/2),q))^4)+12 (((Subscript[ϑ, 1])(-(π/6),q))^4) (((Subscript[ϑ, 1])'(0,(q^2)))^2)=0.

(((Subscript[ϑ, 1])(-(π/3),q))^8)-4 (((Subscript[ϑ, 1])(-(π/6),q))^2) (((Subscript[ϑ, 1])(-(π/3),q))^6)+3 (((Subscript[ϑ, 1])(-(π/6),q))^8)-3 (((Subscript[ϑ, 1])(-(π/2),q))^6) (((Subscript[ϑ, 1])(-(π/6),q))^2)+72 (((Subscript[ϑ, 1])(-(π/6),(q^2)))^2) (((Subscript[ϑ, 1]^′)(-(π/6),(q^2)))^2)=0.

Write a program that, based on the series expansion of the theta functions in q, searches for (and finds) such identities.

θ = EllipticTheta;
Table[Series[θ[j, z, q], {q, 0, 20}], {j, 4}]

For brevity of notation, we will use ϑ for the theta functions. The two functions toTypeset and fromTypeset convert from θ (used for inputting the theta functions) to ϑ (used for formatting the resulting identities).

toTypeset[expr_] := expr //. {θ[a_, z_, q_] :> Subscript[ϑ, a][z, q], 
                              EllipticThetaPrime[a_, z_, q_] :> 
                                                 Derivative[1][Subscript[ϑ, a]][z, q]};
                                                 
fromTypeset[expr_] := expr //. {Subscript[ϑ, a_][z_, q_] :> θ[a, z, q], 
                                Derivative[d_][Subscript[ϑ, a_]][z_, q_] :> 
                                               Derivative[0, 1, 0][θ][a, z, q]};

Generalizing the identities given in the exercise, we will look for identities of the form

(∏_(j=1)^o)(Subscript[c, j])((∏_(k=1)^(Subscript[n, j]))(q^(Subscript[α, j,k]))(((Subscript[ϑ, (Subscript[a, j,k])]^(((Subscript[d, j,k]))))((Subscript[z, j,k]),(q^(Subscript[β, j,k]))))^(Subscript[γ, j,k])))=0

with rational (Subscript[c, k]), (Subscript[z, k,j]) fractions of π, and positive integers (Subscript[α, k,j]), (Subscript[β, k,j]) , (Subscript[γ, k,j]), the order of the derivatives is (Subscript[d, j,k])∈{0,1}, and (Subscript[a, j,k])∈{1,2,3,4}. We restrict the search to identities containing the theta functions and their first derivatives with respect to z. Derivatives with respect to q can be expressed through derivatives with respect to z by using the differential equation 4 q ∂(Subscript[ϑ, a])(z,q)/∂q=-(∂^2)(Subscript[ϑ, a])(z,q)/∂(z^2). And higher derivatives with respect to z can be reduced to the theta functions and their first derivatives (and the single z-independent term (Subscript[ϑ, 4]^″)(0,q)) [✶].

To find such identities programmatically, we form series in q with integer coefficients from the theta functions and their derivatives and then try to form linear combinations of products of such series with unknown coefficients (Subscript[c, j]). Taking the coefficients of the powers of q results in a linear system in the coefficients. Taking enough such linear equations, we will obtain overdetermined systems, and we will search for nontrivial null spaces. The null spaces will then give the desired identities.

Because we want series in q with integer coefficients, but for z=p/q× π we might obtain overall radical factors for the series, we define a function extractCommonRadicalFactor that extracts common radical factors from a series. (We could, of course, in principle use any exact coefficients. But (modular) integer arithmetic is much faster than calculations with arbitrary algebraic numbers.)

extractCommonRadicalFactor[f_, q_] :=
Module[{ser = Series[f, {q, 0, 40}], fac, rad},
       fac = Factor[PowerExpand[Normal[ser] /. q -> 𝓆^4]];
       rad = fac /. (r:Power[_?NumericQ, _Rational]) rest__ :> {r};
       f/Which[Head[rad] =!= List || rad === {}, 1,  True, rad[[1]]]]

We will keep our series in the form {ϑSymbolic,explicitSeriesInq}. For a given expression containing theta functions, the command makeSeriesPair forms such a pair.

makeSeriesPair[f_] := {toTypeset[#], Series[#, {q, 0, 30}]}& @ 
                                     extractCommonRadicalFactor[f, q]

Now we must decide which fractions of π to take for the z-values. The function makeSimpleSeriesTerms form the theta pairs {(Subscript[ϑ, a]^((d)))(z,q),…}, a=1,2,3,4, d=0,1 for a given z.

makeSimpleSeriesTerms[piFractions_] := 
 DeleteCases[#, {_, _?(Normal[#] === 0&)}]& @ 
     Flatten[Table[makeSeriesPair[Derivative[0, d, 0][θ][a, #, q]], 
                   {a, 4}, {d, 0, 1}]& /@ piFractions, 2];

For z=π/8, we cannot extract a simple radical. We form linear combinations of two theta functions instead.

T1 = Join[makeSimpleSeriesTerms[{0, Pi/2, -Pi/2, Pi/4, -Pi/4}],
          makeSeriesPair /@ 
        Join[{θ[3, Pi/8, q] + θ[4, Pi/8, q], θ[3, Pi/8, q] - θ[4, Pi/8, q]},
             {θ[3, Pi/8, q] + θ[4, Pi/8, q], θ[3, Pi/8, q] - θ[4, Pi/8, q]} /.
                                                θ -> Derivative[0, 1, 0][θ]]];
T1 // TableForm // TraditionalForm

We will not only consider (Subscript[ϑ, a]^((d)))(z,q), but the more general term (((Subscript[ϑ, a]^((d)))(z,(q^β)))^γ). So, for the current choices of z, we obtain 96 theta series.

makeBaseTerms[T1_] :=  Flatten[Table[Join[First /@ T1, {} (* First /@ T2 *)] /.
                                            q -> q^(2^e), {e, 0, 2}]];
                                            
ϑBaseTerms = makeBaseTerms[T1];
Length[ϑBaseTerms]

We now calculate the series of the functions ϑBaseTerms to order O((q^250)).

o = 250;

makeϑBaseTermsSeries[ϑBaseTerms_] := DeleteCases[
 {toTypeset[#], ExpandAll @ Series[fromTypeset[#], 
                                   {q, 0, o}]}& /@ ϑBaseTerms,
                               {_, _?((Normal[#] === 0)&)}];

ϑBaseTermsSeries = makeϑBaseTermsSeries[ϑBaseTerms];

Removing all the series that are simple multiples of each other leads to 33 different series. The function identifyIntegerMultiples identities series that are multiples of each other.

(* pull out a common integer from a series *)
pullOutIntegerFactor[{f_, ser_}] := 
Module[{g = Abs[FactorSquareFree[Normal[ser]] /. _Plus :> 1 /. q -> 1]}, 
       If[g == 0, {f, g}, {f/g, ser/g}]]
       
(* make pairs of plus/minus the series *)
sortedPlusMinusPair[{ϑ_, ser_}] := 
Module[{p = {ser, -ser}, pS = Sort[{ser, -ser}]},
       If[p === pS, {pS, {ϑ, -ϑ}}, {pS, {-ϑ, ϑ}}]]

identifyIntegerMultiples[listOfSeries_] := 
Module[{listOfSeries1, listOfSeriesPairs, equivalenceSeries}, 
 (* pull out common integers *)
 listOfSeries1 = pullOutIntegerFactor /@ listOfSeries;
 (* form ±1×series *)
 listOfSeriesPairs = sortedPlusMinusPair /@ listOfSeries1;
 (* identify equivalent series *)
 equivalenceSeries = Split[Sort[listOfSeriesPairs], #1[[1]] === #2[[1]]&];
 (* return equivalence classes in form {ϑ, series} *)
 Apply[{#2[[1]], #1[[1]]}&, equivalenceSeries, {2}]]

equivalentSeriesList = identifyIntegerMultiples[ϑBaseTermsSeries];
ϑBaseTermsFinal = First /@ equivalentSeriesList;
Length[ϑBaseTermsFinal]

Here are the form of the resulting 33 series shown.

Reverse[First /@ ϑBaseTermsFinal] // TraditionalForm

And here are the tuples of theta functions, representing identical series. (For the further calculations, we use only one representative version of these series. By substituting the equivalent forms, we could generate more identities.)

(* pull out an overall factor from the identities if possible *)
overallFactor[l_] := Factor[Plus @@ l] /. _Plus?(Not[NumericQ[#]]&) :> 1

(sIds = Equal @@@ ((#/overallFactor[#])& /@ 
   Map[First, Select[equivalentSeriesList, (Length[#] > 1)&], {2}])) //
                                                        TraditionalForm

A quick numerical high-precision check confirms the last result.

Off[N::meprec];
Union[N[Normal[fromTypeset[sIds]] /. q -> 1/2, 22]]
On[N::meprec];

Some of the series in ϑBaseTermsFinal are powers series in (q^(1/4)) and (q^(1/2)) rather than in q. By multiplying them with each other, we obtain power series in q. The function makeϑBaseProductsPre forms all possible such products.

makeϑBaseProductsPre[ϑBaseTermsFinal_] := 
Module[{fBases2, fBases4, intBases, λ2, λ4},
 (* series in half-integer and in quarter-integer powers of q *)                                    
 {fBases2, fBases4} = 
  Function[q, Cases[ϑBaseTermsFinal, _?(MemberQ[Normal[#[[2]]], 
                    Power[_, Rational[_, q]], Infinity]&)]] /@ {2, 4};                                    
                                     
 (* series in integer powers of q *)
 intBases = Complement[ϑBaseTermsFinal, Join[fBases2, fBases4]];
 (* form all eries in integer powers of q by multiplication *)
 {λ2, λ4} = Length /@ {fBases2, fBases4};
 Join[intBases, 
      (* 1/2 + 1/2 = 1 *)
      Union @ Flatten[Table[((fBases2[[i]]^#1) (fBases2[[j]]^#2))& @@@ 
                                  {{1, 1}}, {i, λ2}, {j, λ2}], 2],
      (* 1/2 + (1/4 + 1/4) = (1/4 + 1/4) + 1/2 = 1 *)
      Union @ Flatten[Table[((fBases2[[i]]^#1) (fBases4[[j]]^#2))& @@@ {{1, 2}},
                            {i, λ2}, {j, λ4}], 2],  
      (* (1/4 + 1/4) + (1/4 + 1/4) 
         = (1/4) + (1/4 + 1/4 + 1/4) = (1/4 + 1/4 + 1/4) + (1/4) = 1 *)
      Union @ Flatten[Table[((fBases4[[i]]^#1) (fBases4[[j]]^#2))& @@@ 
                             {{2, 2}, {1, 3}, {3, 1}},  {i, λ4}, {j, λ4}], 2]]];

By eliminating again doubles, and series that are just negtives of each other from the 73 series, we are left with 72 different series.

ϑBaseProductsPre = makeϑBaseProductsPre[ϑBaseTermsFinal];
ϑBaseProductsPre // Length

equivalentProductSeriesList = identifyIntegerMultiples[ϑBaseProductsPre];
ϑBaseProducts = First /@ equivalentProductSeriesList;
Λ = Length[ϑBaseProducts]

Here are the two products that are equivalent.

Equal @@@ ((#/overallFactor[#])& /@ Take[Map[First, 
   Select[equivalentProductSeriesList, Length[#] > 1&], {2}], All]) // 
                                                       TraditionalForm

And these are the 72 products that we will use as building blocks in forming linear combinations of products and powers to find identities.

(* truncate a series in such a way to display n nonvanishing terms *)
makeNTermSeries[HoldPattern[sd:SeriesData[x_, x0_, cs_, min_, max_, den_]], 
                n_Integer?Positive] := 
Module[{pos = Position[cs, _?(# !=0&), {1}, n], λ}, 
       λ = Max[pos];
       If[Length[pos] < n, sd,
          SeriesData[x, x0, Take[cs, λ], min, min + λ, den]]]

{#1, makeNTermSeries[#2, 4]}& @@@ ϑBaseProducts // TableForm // 
                                                TraditionalForm

Now, we have the theta monomials to use on the ansatz for the linear combinations to form the identities we are looking for. We use all basic monomials with equal probability and raise them to randomly chosen small powers. (Because sometimes the same monomial will be chosen multiple times, this does not limit the maximal power in which a theta function can occur.)

𝕣 := Random[Integer, {1, Λ}];
e𝕣 :=  Random[Integer, {0, 2}];
o𝕣 :=  Random[Integer, {0, 6}];

The function findIdentityCandidates carries out the above-described procedure to find potential nontrivial identities. The argument n indicates the number of theta monomials to be used in the initial ansatz. To speed up the null space calculation, we do not carry it out over rational numbers or integers (which would potentially lead to huge integers and so would be slow), but rather form only a modular null space.

findIdentityCandidates[n_] := 
Module[{prime = Prime[10^6], ansatzFunctions, λ, cList, qPoly, qSeries,
        coeffs, ℓ, eqs, ns, potentialIdentities}, 
(* carry out identity search until a candidate is found *)
While[
 (* form products of powers of theta functions *)
 ansatzFunctions = DeleteCases[#, 1]& @ 
         Union[Table[(* q^e *) Product[ϑS[𝕣]^e𝕣, {o𝕣}], {n}]];
 (* number of different theta monomials *)
 λ = Length[ansatzFunctions];
 (* form linear combination with unknown coefficients *)
 cList = Table[c[k], {k, λ}];
 (* form resulting series *)
 qPoly = cList.ansatzFunctions /. ϑS[k_] :> ϑBaseProducts[[k, 2]];
 qSeries = qPoly + O[q]^o; 
 (* extract linear equations in the c's from q-powers *)
 coeffs = DeleteCases[qSeries[[3]], 0, {1}];
 ℓ = Length[coeffs];
 (* form coefficient matrix with more rows than columns *)
 eqs = Function[sum, Coefficient[sum, #]& /@ cList] /@ 
                                 Take[coeffs, Min[3 λ, ℓ]];
 (* find modular null space *)
 ns = NullSpace[eqs, Modulus -> prime];
 (* lift modular null space *)
 res = (Mod[#, prime, -Floor[prime/2]].ansatzFunctions)& /@ ns;  
 (* the potential identities; ignore multiplicative factors *)
 potentialIdentities = Cases[Factor[res /. ϑS[k_] :> ϑBaseProducts[[k, 1]]], 
                             _Plus, Infinity];
 potentialIdentities === {}];
 (* return found identities *) potentialIdentities]

Because of the use of modular arithmetic in the last function, there is no guarantee that the result holds for all moduli, we add an additional high-precision numeric test to check found identities. The function findIdentity finds a numerically verified identity.

finalqTestQ[id_] := 
Module[{tf},
       Off[N::meprec]; 
       tf = N[fromTypeset[id] /. q -> 1/2, 1000] == 0;
       Off[N::meprec]; 
       tf]

findIdentity[n_] := 
Module[{ids}, 
       While[(ids = Select[findIdentityCandidates[n], finalqTestQ]) == {}];
       ids]

We are now in a position to carry out the identity search. To avoid that we get multiple copies of the same identity, we define a function findIdentities that tries to find o identities. It eliminates all multiple found identities.

findIdentities[o_, n_] :=
Module[{bag = {}, newIds, newIds1, newIds2},
 While[Length[bag] < o,
       newIds = findIdentity[n];
       (* were the newly found identities already encountered *)
       newIds1 = Union[Cases[Factor[newIds], _Plus, Infinity]];
       newIds2 = Select[newIds1, FreeQ[bag, #]&]; 
       (* update identity collection *)
       bag = Join[bag, newIds2]];
 (* return found identities *) bag]

And here are 25 identities involving some of the above-constructed 33 theta terms.

SeedRandom[2];
findIdentities[25, 300]

By changing the values of z as fractions of π to another value uniformly, we can extract the identities that are valid for any z. (To find such identities, we could have started with series expansions in q and z and form systems of equations in powers of q and z.)

Select[fromTypeset[%], (N[# /. Pi -> E /. q -> 1/2, 22] == 
                        N[# /. Pi -> 1/Pi /. q -> 1/2, 22] == 0)&]

We can now easily carry out a search for theta identities with different values of z. Here is another list of theta function series.

T1 = makeSimpleSeriesTerms[{0, Pi/2, -Pi/2, Pi/3, -Pi/3, Pi/6, -Pi/6}];

This time, we have 36 theta monomials.

ϑBaseTerms = makeBaseTerms[T1];
ϑBaseTermsSeries = makeϑBaseTermsSeries[ϑBaseTerms];
equivalentSeriesList = identifyIntegerMultiples[ϑBaseTermsSeries];
ϑBaseTermsFinal = First /@ equivalentSeriesList;
Length[ϑBaseTermsFinal]

And for the linear combinations, we can draw from 132 theta monomials and no equivalent series among the ϑBaseProductsPre.

ϑBaseProductsPre = makeϑBaseProductsPre[ϑBaseTermsFinal];
equivalentProductSeriesList = identifyIntegerMultiples[ϑBaseProductsPre];
ϑBaseProducts = First /@ equivalentProductSeriesList;
Λ = Length[ϑBaseProducts]

And here are again 20 identities among them.

SeedRandom[2];
findIdentities[20, 360]

The interested reader could now continue with further theta series. Here are a few further possibilities for the series to be used in ϑBaseTerms.

makeSeriesPair /@ Join[
   { θ[1, Pi/12, q] + θ[2, Pi/12, q], θ[1, Pi/12, q] - θ[2, Pi/12, q],
     θ[3, Pi/12, q] + θ[4, Pi/12, q], θ[3, Pi/12, q] - θ[4, Pi/12, q],
     θ[2, 4 Pi/5, q] - θ[2, 2 Pi/5, q],
    -θ[2,  Pi/5, q] + θ[2, 2 Pi/5, q] - θ[2, 3 Pi/5, q] + θ[2, 4 Pi/5, q],
     θ[3, 3 Pi/5, q] + θ[3, 4 Pi/5, q],
    -θ[3, 3 Pi/10, q] - θ[3, 7 Pi/10, q] + θ[3, 9 Pi/10, q] + 3 θ[4, Pi/5, q],
    -θ[3,   Pi/10, q] - θ[3, 3 Pi/10, q] - θ[3, 7 Pi/10, q] + 3 θ[4, Pi/5, q]},
   { θp[1, Pi/12, q] + θp[2, Pi/12, q], θp[1, Pi/12, q] - θp[2, Pi/12, q],
     θp[3, Pi/12, q] + θp[4, Pi/12, q], θp[3, Pi/12, q] - θp[4, Pi/12, q]} /.
                                                  θp -> Derivative[0, 1, 0][θ]]

These are some of the series that result from the use of theta functions at values that are multiples of π/5 and π/10.

-(Subscript[ϑ, 2])((π/5),(q^4))+(Subscript[ϑ, 2])(((2 π)/5),(q^4))-(Subscript[ϑ, 2])(((3 π)/5),(q^4))+(Subscript[ϑ, 2])(((4 π)/5),(q^4))+2 (-(Subscript[ϑ, 3])(((3 π)/5),q)+(Subscript[ϑ, 3])(((3 π)/5),(q^4))-(Subscript[ϑ, 3])(((4 π)/5),q)+(Subscript[ϑ, 3])(((4 π)/5),(q^4)))=0

(Subscript[ϑ, 2])((π/5),(q^4))-(Subscript[ϑ, 2])(((2 π)/5),(q^4))+(Subscript[ϑ, 2])(((3 π)/5),(q^4))-(Subscript[ϑ, 2])(((4 π)/5),(q^4))+2 ((Subscript[ϑ, 3])(((3 π)/5),(q^4))+(Subscript[ϑ, 3])(((4 π)/5),(q^4))+(Subscript[ϑ, 3])(((3 π)/10),q)+(Subscript[ϑ, 3])(((7 π)/10),q)-(Subscript[ϑ, 3])(((9 π)/10),q)-3 (Subscript[ϑ, 4])((π/5),q))=0

-(Subscript[ϑ, 3])(((3 π)/5),q)+2 (Subscript[ϑ, 3])(((3 π)/5),(q^4))-(Subscript[ϑ, 3])(((4 π)/5),q)+2 (Subscript[ϑ, 3])(((4 π)/5),(q^4))+(Subscript[ϑ, 3])(((3 π)/10),q)+(Subscript[ϑ, 3])(((7 π)/10),q)-(Subscript[ϑ, 3])(((9 π)/10),q)-3 (Subscript[ϑ, 4])((π/5),q)=0

-(Subscript[ϑ, 2])((π/5),(q^4))+(Subscript[ϑ, 2])(((2 π)/5),(q^4))-(Subscript[ϑ, 2])(((3 π)/5),(q^4))+(Subscript[ϑ, 2])(((4 π)/5),(q^4))-(Subscript[ϑ, 3])(((3 π)/5),q)-(Subscript[ϑ, 3])(((4 π)/5),q)-(Subscript[ϑ, 3])(((3 π)/10),q)-(Subscript[ϑ, 3])(((7 π)/10),q)+(Subscript[ϑ, 3])(((9 π)/10),q)+3 (Subscript[ϑ, 4])((π/5),q)=0

((((Subscript[ϑ, 3])(((3 π)/5),q)+(Subscript[ϑ, 3])(((4 π)/5),q)))^2)-((((Subscript[ϑ, 3])(((3 π)/10),q)+(Subscript[ϑ, 3])(((7 π)/10),q)-(Subscript[ϑ, 3])(((9 π)/10),q)-3 (Subscript[ϑ, 4])((π/5),q)))^2)+2 ((Subscript[ϑ, 2])((π/5),(q^4))-(Subscript[ϑ, 2])(((2 π)/5),(q^4))+(Subscript[ϑ, 2])(((3 π)/5),(q^4))-(Subscript[ϑ, 2])(((4 π)/5),(q^4))) ((Subscript[ϑ, 3])(((3 π)/5),(q^4))+(Subscript[ϑ, 3])(((4 π)/5),(q^4)))=0

2 (((((Subscript[ϑ, 3])(((3 π)/5),q)+(Subscript[ϑ, 3])(((4 π)/5),q)))^2)-2 ((((Subscript[ϑ, 3])(((3 π)/5),(q^4))+(Subscript[ϑ, 3])(((4 π)/5),(q^4))))^2)+((((Subscript[ϑ, 3])(((3 π)/10),q)+(Subscript[ϑ, 3])(((7 π)/10),q)-(Subscript[ϑ, 3])(((9 π)/10),q)-3 (Subscript[ϑ, 4])((π/5),q)))^2))-((((Subscript[ϑ, 2])((π/5),(q^4))-(Subscript[ϑ, 2])(((2 π)/5),(q^4))+(Subscript[ϑ, 2])(((3 π)/5),(q^4))-(Subscript[ϑ, 2])(((4 π)/5),(q^4))))^2)=0

To search for identities with other multiples of π/d, we need linear combinations of (Subscript[ϑ, a])(p π/d,q) such that the series expansion of these linear combinations has integer coefficients. The following function findϑLinearCombinations tries to find o linear combinations of thsi form for multiples of π/den.

findϑLinearCombinations[o_, den_] := 
Module[{(* general series expansions at q == 0 *)
        ser = Normal[Table[Series[EllipticTheta[j, z, q], 
                                  {q, 0, 25}], {j, 4}]],
        sum12, sumS12, sum34, sumS34, vars, 𝕣𝕔, 𝕣𝕔𝕤, 𝕤, 
        n𝕤, n𝕤r, bag = {}},
 (* symbolic linear combinations of the two pairs with equal q-powers *)
 {sumS12, sumS34} = 
   (Sum[c[α, j] Subscript[ϑ, #1][j Pi/(2 den), q], {j, 0, 2 den - 1}] + 
    Sum[c[β, j] Subscript[ϑ, #2][j Pi/(2 den), q], {j, 0, 2 den - 1}])& @@@
                                                           {{1, 2}, {3, 4}};
 (* explicit linear combinations of the two pairs with equal q-powers *)
 {sum12, sum34} = {sumS12, sumS34} /. 
                Subscript[ϑ, j_] :> Function @@ {{z, q}, Evaluate[ser[[j]]]};      
 (* the undetermined coefficients *)
 vars = Cases[sumS12, _c, Infinity];
 (* values for the 4 den undetermined coefficients 
    c[α, j] and c[β, j], j = 0, 1, …, 2 den - 1 *)
 𝕣𝕔 := Rule @@@ Transpose[{vars, Table[Random[Integer, {-1, 1}], 
                                       {Length[vars]}]}];
 While[Length[bag] < o, (* try linear combination *) 𝕣𝕔𝕤 = 𝕣𝕔;       
       (* check if try combination results in q-series 
       with integer coefficients *)     
       (* series in ϑ_1 and ϑ_2 *)                     
       𝕤 = sum12 //. 𝕣𝕔𝕤; n𝕤 = N[𝕤, 22]; n𝕤r = n𝕤 /. x_Real :> Round[x];
       If[(n𝕤 - n𝕤r /. q -> 1) == 0, AppendTo[bag, {sumS12 /. 𝕣𝕔𝕤, n𝕤r}]];
       (* series in ϑ_3 and ϑ_4 *) 
       𝕤 = sum34 //. 𝕣𝕔𝕤; n𝕤 = N[𝕤, 22]; n𝕤r = n𝕤 //. x_Real :> Round[x];
       If[(n𝕤 - n𝕤r  /. q -> 1) == 0, AppendTo[bag, {sumS34 /. 𝕣𝕔𝕤, n𝕤r}]]];
 (* return found linear combinations *) bag]

Here are some linear combinations of (Subscript[ϑ, a])(p π/7,q) that have integer q-series coefficients.

SeedRandom[123];
findϑLinearCombinations[6, 7]