Rod-packing arrangements of invariant tori in solenoidal vector fields with cubic symmetries

The arrangements of invariant tori that resemble rod packings with cubic symmetries are considered in three-dimensional solenoidal vector fields. To find them systematically, vector fields whose components are represented in the form of multiple Fourier series with finite terms are classified using magnetic groups. The maximal magnetic group compatible with each arrangement is specified on the assumption that the cores of the nested invariant tori are straight and located on the lines corresponding to the central axes of the rods packed. Desired rod-packing arrangements are demonstrated by selecting vector fields whose magnetic groups are the maximal ones and by drawing their integral curves that twine around invariant tori. In the demonstration of chiral arrangements, Beltrami flows (or force-free fields in plasma physics), which have the strongest chirality of all solenoidal vector fields satisfying the same vector Helmholtz equation, are used. As by-products, several chain-like arrangements of closed invariant tori were found. One of the chains consists of knotted invariant tori. In all vector fields (chiral or achiral) selected for the demonstration, the volume percentages of ordered regions formed by invariant tori in a unit cell were roughly measured with the aid of a supervised machine learning technique.


Introduction
The problem of sphere packings has a long history since Kepler [1] mentioned it in 1611 after becoming aware of the work of T. Harriot (see [2] for the history and the proof of the so-called Kepler conjecture on dense sphere packings). By contrast, the problem of rod (or cylinder) packings is much newer. The pioneering work of O'Keeffe and Andersson [3] in inorganic chemistry was published in 1977. Recently, rod packings have attracted the interest of chemists in the design of materials such as metal-organic frameworks or covalent organic frameworks [4,5]. The rod-packing problem was generalised by assuming that rods are not necessarily straight but may be zigzagged, curved or tangled [5,6].    Various rod packings exist even when only considering those with cubic symmetries [7]. However, if we restrict ourselves to invariant cubic rod packings, then the number of packings is six (or nine if a pair of enantiomorphs is counted as two). Here, an invariant cubic rod packing means a packing in which the axes of rods coincide with the rotation or screw rotation axes (with or without inversion centres) specified in a cubic space group. The six (or nine) packings were named Π , Σ , Ω (or + Π , − Π , + Σ , − Σ , + Ω , − Ω ), Π * , Σ * and Γ in [8] (see also [5,Fig. 16]), and their rods are arranged as shown in Figs. 1, 2, 3, 4, 5, and 6. Here, the signs ± distinguish enantiomorphs, and the sign * indicates the union of the enantiomorphs.
O'Keeffe [9] referred to the Π and Π * packings as -Mn and -W packings, respectively, because they are convenient for describing the crystal structures of -manganese and -tungsten. For the same reason, the Σ , Σ * and Γ packings (except Ω ) were referred to as the SrSi 2 , -Si and garnet packings, respectively. The Π and Π * packings also arise in blue phases I and II, respectively, of cholesteric liquid crystals. They were theoretically proposed in the 1980s [10,11] and experimentally observed approximately 30 years later [12]. The rods in the blue phases are called double-twist cylinders. Each double-twist cylinder consists of a nested set of coaxial cylindrical surfaces around which liquid crystalline molecules wrap helically (see e.g. [13,14]).
Intriguingly, a figure similar to that in Fig. 1 was presented in a fluid mechanics paper [15], which investigated the structure of integral curves (or streamlines) of the Arnold-Beltrami-Childress (ABC) flow with arbitrary real constants A, B and C. The figure [15,Fig. 4] sketched the arrangement of tubes, consisting of nested tubelike surfaces helically wound around by integral curves of with ABC ≠ 0 . In studies of dynamical systems, such a tubelike surface is referred to as an invariant torus [16,17]. When |A| = |B| = |C| ≠ 0 , the ABC flow acquires cubic symmetry. Figure 7a shows a family of its invariant tori, which were called principal vortices in [15]. Apart from their undulations, they are arranged similarly to the rods in the − Π packing in Fig. 1 Fig. 7a is depicted by numerically calculating the integral of :

. Each invariant torus in
which is the solution to the dynamical system = 0 , and plotting the curve of (t) modulo the unit cell [ ∕2, 5 ∕2) 3 . The flow with |A| = |B| = |C| ≠ 0 also has invariant tori of another family, which were called secondary vortices in [15]. Two of these are shown in Fig. 7b, in which an eight times larger cell is used to clearly show that the two tori form a double helix 1 along a rod in the + Σ packing. Hereafter, for brevity, we refer to the arrangement of invariant tori along the rods in the + Π packing, etc. as the + Π arrangement, etc.
In general, invariant tori in three dimensions arise in the integration of solenoidal (i.e. volume-preserving) vector fields. Indeed, satisfies div = 0 . If an initial point is chosen outside the invariant tori (and moreover, it is neither a zero of the vector field nor located on a heteroclinic trajectory connecting two zeros), then the integral curve chaotically wanders without belonging to any invariant torus, as shown for in Fig. 8. 2 In other words, space is divided into an ordered part with invariant tori and a chaotic part with chaotic integral curves. The former and the latter are sometimes called "islands" and a "sea", respectively, by their appearances on the Poincaré section (e.g. [18, p. 427]). For example, with ABC = 0 is analytically integrable (see [15]), and there exists no "sea". By contrast, in with |A| = |B| = |C| ≠ 0 , the volume percentages of the above mentioned − Π -and + Σ -arranged ordered regions and the chaotic region per unit cell are estimated to be 87.5%, 5.6% and 6.9%, respectively, using the method explained in Sect. 2. In every vector field investigated in this study, a nonzero percentage of a unit cell is chaotic. We should note that chaos in solenoidal vector fields is different from chaos in Fig. 7 Invariant tori in the ABC flow with A = B = C = 1. a Arrangement similar to the − Π rod packing. The six integral curves that start from 0 = (1.9 , 0.8 , 0.5 ) , (0.9 , 1.8 , 1.5 ) and their respective cyclic permutations are plotted modulo the unit cell [ ∕2, 5 ∕2) 3 . b Arrangement along + Σ-packed rods. Only the two integral curves that start from 0 = (0.9 , 0.7 , 0.5 ) and (1.9 , 1.7 , 3.5 ) are plotted modulo the eight times larger cell [ ∕2, 9 ∕2) 3 . A rod in the + Σ packing is plotted in the same way non-solenoidal vector fields such as Lorenz's volume-contracting flow, which is well known for its chaos attractor.
In the late 1980s, a research group in the Soviet Union, e.g. [19,20], investigated nearly two-dimensional solenoidal vector fields that have crystal or quasicrystal structures. Because each of these vector fields is not integrable but nearly integrable, its chaotic integral curves are confined in a thin region, which was called a stochastic web. The cross section of this thin region shows a periodic or quasiperiodic pattern. By contrast, in a fully three-dimensional solenoidal vector field without integrability, like ours, the chaotic region is not usually thin and its complement, the ordered region with invariant tori, is easier to study.
Noticing the resemblance between the Π arrangements of invariant tori in and of double-twist cylinders in blue phase I, we [21] derived three solenoidal fluid flows associated with cubic space groups that had been used in studies of blue phases. The flows were shown to have Γ -, − Ω -and Π * -arranged invariant tori [21,Figs. 4,8 and 11], although the connection with the rod packings was not mentioned. In addition, the first flow has + Σ-arranged invariant tori, as well [21,Fig. 7a]. We [22] extended our study to solenoidal fluid flows with hexagonal symmetries and reported that there exist various arrangements of invariant tori, two of which resemble the Λ and Δ rod packings (see [5,Fig. 16] for these hexagonal packings).
One of the aims of this study is to investigate the type of solenoidal vector field (not restricted to fluid flow) having invariant tori whose arrangement resembles each of the six (or nine) cubic rod packings. For simplicity, we focus on the cases in which the cores of nested invariant tori are straight and located on the same lines as the central axes of the rods. Following the procedure used in [22], we apply magnetic groups to vector fields 3 klm = (u, v, w) such that each of u, v and w is written in the form of multiple Fourier series with at most 48 ( = 3! ⋅ 8 ) terms: Here, a pqri ( i = 1, 2, … , 8 ) are constants, k, l and m are non-negative integers such that (k, l, m) ≠ (0, 0, 0) and (k∕n, l∕n, m∕n) ∉ ℤ 3 for any integer n ≥ 2 , and P(k, l, m) is the set of all permutations of (k, l, m). We will restrict the coefficients a pqri by specifying the magnetic group of klm , and clarify which magnetic group is compatible with each rod-packing arrangement of invariant tori with straight cores.
The other aim of this study is to demonstrate the desired arrangements of invariant tori. Unlike rods in the rod packings, two interlaced sets of nested invariant tori cannot touch because of the uniqueness of the integral curve to a given initial point, and their gap is filled with chaotic integral curves (in the "sea" of chaos). Moreover, even if the cores of invariant tori are straight, the shapes and sizes of the cross sections of off-core invariant tori are nonconstant. Therefore, the volume percentage of invariant tori in each rod-packing arrangement is different from (generally, lower than) the percentage of rods packed. We will show how different they are in each vector field used as an example.
When klm is solenoidal, we write it as klm,S . Because every klm has the x, y and z components in the form (1) and satisfies the vector Helmholtz equation with the solenoidal field klm,S is also a solution to (2). Consequently, the vector field defined as = klm,S exp(i c 0 t) with a constant c 0 > 0 satisfies the vector wave equation as well as the equation div = 0 . This means that can be realised in an electric or a magnetic field by superposing monochromatic electromagnetic waves. Therefore, the study of the arrangements of invariant tori has the potential to be useful in materials science. For example, if we can orientate rod-like molecules along the electric or magnetic force lines and somehow maintain their orientation, then the construction of crystals with artificial double-twist cylinders will be possible. (1) u(x, y, z) = ∑ (p,q,r)∈P(k,l,m) (a pqr1 cos px cos qy cos rz + a pqr2 cos px cos qy sin rz + ⋯ + a pqr8 sin px sin qy sin rz), v(x, y, z) = u(y, z, x), w(x, y, z) = u(z, x, y).
As vector fields with chirality, we define klm,LB and klm,RB by restricting the constants a pqri of klm in (1) so that respectively, are satisfied. These are special cases of klm,S (indeed, div klm,LB = −1 div rot klm,LB = 0 etc.). The ABC flow is an example of 001,LB . The labels LB and RB represent the left-and right-handed Beltrami flows, respectively. In general, a solenoidal vector field that parallels nontrivial rot is called a Beltrami flow (or a force-free field in plasma physics). At the point where the ratio of rot to is positive (resp. negative), locally twists like a left-handed (resp. right-handed) helix. Indeed, the simplest Beltrami flows LB = (sin 0 z, cos 0 z, 0) ( = −1 0 rot LB ) and RB = (cos 0 z, sin 0 z, 0) ( = − −1 0 rot RB ) with a positive constant 0 twist left-and right-handed helically, respectively, at all points. See [23] for details on the local geometry of Beltrami flows.
Because every klm,S is a solution to (2), it also satisfies Using this fact, we can verify that an arbitrary non-Beltrami field klm,S is represented by the sum of and these two fields are left-and right-handed Beltrami flows, respectively (as shown by calculating their curls). See [24], in which the same fact was mentioned in terms of the Fourier series. The decomposition of klm,S into klm,LB and klm,RB is analogous to the decomposition of linearly or elliptically polarised light into leftand right-handed circularly polarised light. For this reason, we can regard klm,LB and klm,RB to have the strongest left-and right-handed chirality, respectively, of all klm,S . A remarkable difference between klm,LB and klm,RB appears in their integral curves, which wind around invariant tori (if existing) helically right-and lefthanded, respectively, with the handedness opposite to that mentioned above (see Fig. 9). As in our previous papers [21,22], we will mainly use klm,RB , rather than klm,LB , for chiral arrangements of invariant tori. For achiral arrangements, we will use klm,S , which are represented with both klm,LB and klm,RB (as it were, analogously to racemic bodies in chemistry). In the achiral arrangements, the type of invariant torus in Fig. 9 and its mirror image coexist.
To integrate vector fields and draw their integral curves as in Fig. 7 or Fig. 8, we used the command "NDSolve" in the software Mathematica. To estimate the volume percentages of ordered or chaotic regions, we applied the machine learning command "Classify" (with the option "Method" designated as "Logis-ticRegression") to a huge quantity of figures showing integral curves whose initial points are equally distributed in a unit cell. After that, we counted the figures of curves in each region by the command "Count", considering their ratio to all figures to approximate the volume percentage of the region. In Sect. 2, details of how to find the rod-packing arrangements of invariant tori are explained after  introducing magnetic groups. In Sect. 3, the Π and Σ arrangements are considered using some magnetic groups of the I4 1 32 family. In Sect. 4, the Σ * arrangement is found in a vector field whose magnetic group belongs to the Ia3d family, which is the superfamily of the I4 1 32 family. In Sect. 5, a magnetic group of the Ia3 family is used for the Π * arrangement. In addition to this achiral Π * arrangement, a chiral Π * arrangement is found using a magnetic group of the I4 1 32 family in Sect. 5, as well.
The I4 1 32 family is also applied to the Γ arrangement in Sect. 6. The Ω arrangement is treated with a magnetic group of the I432 family in Sect. 7. Lastly, in Sect. 8, the results of this study are summarised and discussed.

Magnetic groups
As demonstrated in [22], magnetic groups are convenient to derive vector fields whose sets of integral curves have desired symmetries. In the three-dimensional case, there are 230, 230 and 1191 magnetic groups of the first, second and third types, respectively [25]. Each of the third-type magnetic groups consists of isometry transformations { | } and { | } � (with , ∈ O(3) and , ∈ ℝ 3 ) that map to + and + , and work as symmetry and antisymmetry operations, respectively, of an axial vector field : By contrast, each first-type magnetic group consists of symmetry operations { | } only, and is denoted by the same international symbol with a space group. Attention is required when a first-or third-type magnetic group is applied to a polar vector field. Because the symmetry and antisymmetry operations of a polar vector field are defined by (5) with the factors det and det removed, every element { | } with det = −1 (resp. { | } � with det = −1 ) of the magnetic group works as an antisymmetry (resp. a symmetry) operation of the polar vector field. In this study, all vectors are assumed to be axial so that we can use first-or third-type magnetic groups without such attention. If is equated with − , like a director field in liquid crystals, then its magnetic group is a second-type group, which is the direct product Each of the total 1651 magnetic groups was numbered in the form a.b.c by Litvin [26]. Here, a ( = 1, 2, … , 230 ) is the international number of the associated space group, b ( = 1, 2, 3, … ) classifies the magnetic groups with the same a, and c ( = 1, 2, … , 1651 ) is a serial number. The first-and second-type magnetic groups have b = 1 and 2, respectively, whereas the third-type groups have b ≥ 3 . We denote klm with (1) whose magnetic group number is a.b.c by (a.b) klm (without using c). In the same way, klm,S , klm,LB and klm,RB are written as

001,LB
because its magnetic group I P 4 1 32 is numbered as 214.4.1570 in [26].

Process
The rod-packing arrangements of invariant tori are found according to the following process: 1. Let G be the following space group (or a cubic subgroup 4 of it when returned from step 7): I4 1 (1) klm , using equivalent points to (x, y, z), especially those yielded by the "generators selected" of G in [26]. Because the coordinates of the equivalent points were listed modulo 1 in [26], convert them into those modulo 2 .
klm,RB , obtained in step 4, or a linear combination of vector fields with the same (a.b) and the same label ( S or RB ) but different of In the last case, a linear combination whose have rational ratio is preferable from a physical viewpoint (e.g. in electromagnetic waves mentioned in Sect. 1).

Fig. 9
A sketch of the field twist and the windings of integral curves on invariant tori in a right-handed Beltrami flow klm,RB 6. Specify all constants contained in . If the cores of invariant tori are required to be straight, give appropriate values to the constants so that the orthogonal component of is zero on the line corresponding to the central axis of a rod in the rod packing. 7. Find an initial point 0 such that the integral curve of , i.e. the solution curve (t) of with (0) = 0 , winds around an invariant torus when numerically drawn modulo a unit cell until an appropriately given upper bound t = t max . If a desired invariant torus cannot be found, return to step 6 and respecify the constants, or start again from step 2 with another G . If no desired invariant torus is found for any G of the same a, return to step 1. 8. Draw all the other invariant tori by taking equivalent points to 0 as new initial points in step 7.
We should note that, if (5) is satisfied, then (6) with replaced by { | } and by is the same as of − (if the same (0) is given), both { | } and { | } � preserve the set of all integral curves of . Moreover, { | } also preserves the set. Therefore, the space group of the set is G, which is obtained by changing all antisymmetry operations to symmetry operations (as { | } � to { | } ) in G . This is why the above process works.
Assume that the cores of invariant tori are straight and one of them coincides with the line = + s , where ∈ ℝ 3 is a fixed position vector, s ∈ ℝ is a parameter and is the direction vector defined as Then, we notice that for the Σ, Ω, Σ * or Γ arrangement.
(7) ≠ and ∥ everywhere on the line = + s (see Fig. 9). In other words, if (a.b) klm obtained in step 3 satisfies or then every in step 6 that satisfies | × ( + s )| = 0 for all s also satisfies ⋅ ( + s ) = 0 for some s, meaning (7) is impossible; thus, the line = + s cannot become a core of invariant tori. Using (8) or (9) as a sufficient condition for the nonexistence of invariant tori with straight cores, we can judge whether the selection of G in step 2 is inappropriate or not. Also, using the negation of (8) or (9) as a necessary condition for the existence of the desired invariant tori, we can judge whether G is worth investigating.
In this study, body-centred cubic space groups are used as G. In such a bodycentred cubic case, the volume percentage of an ordered region with invariant tori or of a chaotic region is roughly measured by the following simple method: i. Make 4N 3 figures of integral curves starting from the equally distributed initial Here, n x and n y are integers from 1 to 2N, whereas n z is an integer from 1 to N. The small constants x , y and z are for deviating 0 from the symmetry axes or planes on which zeros or heteroclinic trajectories are located. We set N = 20 , ii. Count the figures with integral curves winding around invariant tori of the same family or the figures with chaotic integral curves. If the number is M, then M∕(4N 3 ) is an approximate value of the volume fraction of the ordered or chaotic region. Because 4N 3 is a huge number when N = 20 , we classified all the 4N 3 figures into several categories according to (dis)orderedness or the shapes of invariant tori by a supervised machine learning technique (and personally corrected misclassifications) before counting.
The reason for taking 4N 3 initial points, not (2N) 3 , in step i is that the unit cube [0, 2 ) 3 has twice the volume of a primitive unit parallelepiped when G is bodycentred cubic.
klm ( + s * ) = 0 with a constant s * for all possible (k, l, m) In the appropriate setting of the origin and the x, y and z axes, the conditions for Here, (−x + , −y, z + ) and (y + ∕2, x + 3 ∕2, −z + 3 ∕2) are the equivalent points to (x, y, z) yielded by Nos. 2 and 13 of the "generators selected" of each magnetic group in [26], whereas (x + , y + , z + ) is yielded by t(1/2, 1/2, 1/2) or t � (1∕2, 1∕2, 1∕2) of the generators. Due to the forms of v and w in (1), the equivalent points yielded by Nos. 3 and 5 of the generators are not necessary. We should note that, if all of k, l and m are odd, then cubic klm has the face-centred structure: klm (x + , y + , z) = klm (x, y, z) , and its magnetic group is not of the I4 1 32 family. Moreover, because (1) has the same form even if k, l and m are permuted, k and m can be assumed to be even and odd, respectively, without generality loss. On this assumption, we introduce the following four cases: Table 1. In the table, the following are referred to: [k ≡ l(≡ 0 or 2) and m ≡ 1 mod 4] or (k ≡ 0, l ≡ 2 and m ≡ 3 mod 4), Here, the constants c i ( i = 1, 2, … , 12 ) are generally nonzero, although some can vanish unless u becomes what will appear in Sect. 4 klm,RB in step 4 of the process are also listed in Table 1 by referring to (18) u = c 1 (cos kx cos ly sin mz ∓ cos kx cos my sin lz) + c 2 (sin kx sin ly cos mz ± sin kx sin my cos lz) + c 3 (cos lx cos my cos kz ∓ sin lx cos ky sin mz) + c 4 (sin lx sin my sin kz ± cos lx sin ky cos mz) + c 5 (cos mx sin ky cos lz ± sin mx sin ly sin kz) + c 6 (sin mx cos ky sin lz ∓ cos mx cos ly cos kz), (20) where and are generally nonzero constants, and is from (3). In the case of (14) with (k, l, m) = (0, 1, 1) , excluded in Table 1, the x component (18) with the upper signs is not of (214.1) klm , but of (230.1) klm , which is a gradient field associated with the "gyroid" as pointed out in [21,Sect. 4.1]. In addition, (21) with (k, l, m) = (0, 0, 1) is also excluded in Table 1 because (214.4) klm,RB becomes trivial. We can easily show that, if k = 0 or k = l or l = m , then each of (214.b) klm,RB ( b = 1, 3, 4, 5 ) is determined uniquely up to constant multiples. We can also show that (22) and (k, l, m) ≠ (0, 0, 1) in (24).

19
(see [28, Fig. 1.5]). Two of the linked knots are depicted in Fig. 10c. The chain of the knots occupies approximately 11.1% of a unit cell.
We next show that the − Π arrangement of invariant tori with straight cores is realised (not only with

041,RB
, although the ratio of their wavelengths is not rational. They are obtained from Table 1 u(x, y, z) =2 cos 2x(cos 2y cos z + sin y cos 2z) + sin 2x(sin 2y sin z + cos y sin 2z) − 2(cos x cos 2y sin 2z − sin x sin 2y cos 2z) u(x, y, z) = √ 5(cos 2y cos z + sin y cos 2z) + cos 2x(sin y + cos z) − 2(cos x sin 2z − sin x sin 2y), with constant multiples appropriately set. From now on, vector fields unique up to constant multiples will be introduced with appropriate values given to the constants, which are independent of the shapes of invariant tori. Each of (28) and (29) has undulating invariant tori along the 4 3 screw axes. When the linear combination = (214.4) 021,RB + (214.4) 041,RB is proven to satisfy (7) on the 4 3 screw axes (e.g. (3 ∕2, 0, z) with any z) and expected to have − Π-arranged invariant tori with straight cores. Indeed, these tori are observed as shown in Fig. 11. The volume percentage of the ordered region formed by them is estimated at 5.3%. As explained by Fig. 9 in Sect. 1, integral curves in a right-handed Beltrami flow wind around invariant tori left-handed helically. This left-handedly helical winding is confirmed in Fig. 11 as well as in Fig. 10a (cf. Fig. 7a).
From Table 1, we obtain This has thick undulating invariant tori along the 3 1 screw axes, as shown in [21, Fig. 7a]. To cancel its orthogonal component on the 3 1 or 3 2 screw axes, we make use of the following 1 5 -wavelength field:    (29) and (30). One of the initial points is (10) with (n x , n y , n z ) = (29, 2, 1)

Note that (214.b)
klm,RB is not unique even up to constant multiples when k, l and m are nonzero and different from one another. The constants and for which the sum of (32) and (33) (7) is not satisfied. Then, we additionally use (214. 3) 011,RB (3x, 3y, 3z) as a field of period 2 in x, y and z. Because it satisfies not only the latter of (4) with = √ 3 2 + 3 2 but also (11)-(13) for b = 3 , we write it as (214. 3) 033,RB (x, y, z) , i.e.
which is not given by Table 1. The reason for introducing (35) (a 1 3 -wavelength field of (214.3)

011,RB
) is that its orthogonal component is zero and its parallel component is nonzero everywhere on the 3 1 (and also the 3 2 ) screw axes. Indeed, Therefore, can be selected without changing (34) so that (7)  033,RB and + Σ-arranged invariant tori with straight cores are found. In fact, such tori exist as shown in Fig. 12a, where = −9 . The volume percentage of the ordered region formed by them is estimated at 1.4% . This is lower than the percentage of rods in the + Σ or − Σ packing ( √ 3 ∕72 = 7.55 … % , see [7,9] For these and , the orthogonal component of the vector field vanishes on the 3 2 screw axes independently of . When we set = 3 , thin − Σ-arranged invariant tori, which occupy approximately 0.7% of a unit cell, are observed as in Fig. 12b ( b = 1, 3, 4, 5 ) is to satisfy (11) and the following three equations: The four equations (11) and (37)-(39) correspond to Nos. 2, 13, 25 and t(1/2, 1/2, 1/2), respectively, of the "generators selected" of each magnetic group in [26]. For the same reason as in Sect. 3.1, we assume that k and m are even and (37)    (32), (33) and (35). a The + Σ arrangement when and are given by (34) and = −9 . One of the initial points is (10) with (n x , n y , n z ) = (13, 28, 1) . b The − Σ arrangement when and are given by (36) and = 3 . One of the initial points is (10) with (n x , n y , n z ) = (28,14,1) odd, respectively. Moreover, (39) implies that k + l + m is even; thus, we restrict ourselves to the cases (14) and (15). Comparing (37) and (39) with (12) and (13) There are two places where the case (k, l, m) = (0, 1, 1) is excluded in Table 2, for the vector fields become trivial.
We discuss the case in which G is a space group of the first (achiral) category in Sect. 5.2 after studying (206.b) klm,S in Sect. 5.1. The case in which G belongs to the second (chiral) category is treated in Sect. 5.3. These conditions are yielded by the "generators selected" of each magnetic group in [26]. As in the previous sections, we assume that k is even and m is odd. The general forms of the x components u of (206.b) klm are tabulated in Table 3, in which the following are referred to: (52) u = C 1 sin kx sin ly cos mz + C 2 sin kx sin my cos lz + C 3 cos lx cos my cos kz + C 4 sin lx cos ky sin mz + C 5 sin mx cos ky sin lz + C 6 cos mx cos ly cos kz,

Achiral 5 * arrangement
As mentioned at the beginning of Sect. 5, we start the process taking Ia3d as G. According to the columns " -W (×2) " in [9, Table 1], the central axis of a rod in the Ia3d-symmetric Π * packing coincides with the line (0, ∕2, z) or ( ∕2, 0, z) in our coordinates. In this case, the eight unit cells in Fig. 4 becomes one unit cell, as indicated by " (×2) " in [9, Table 1].
In addition, the line ( ∕2, 0, z) in (230.b) klm,S cannot become a core of invariant tori, either. Indeed, at the point ( ∕2, 0, 0) on the line, w( ∕2, 0, 0) ( = u(0, ∕2, 0) ) is equal to zero in (40) and (41), namely (8) holds. Thus, the Π * arrangement of straight invariant tori cannot be expected when Ia3d is chosen as G. The maximal non-isomorphic subgroups of Ia3d are I43d , I4 1 32 and Ia3 . The first one is discarded in the same way as Ia3d . The second one will be taken as G in Sect. 5.3.
When the third one, Ia3 , is taken as G, none of the resulting fields (206.b) klm,S (see Table 3) has a core of invariant tori on the line (0, ∕2, z) . This is because (61) remains valid even if klm,S has a core of invariant tori on the line ( ∕2, 0, z) . For these reasons, using (206.4) klm,S , we search for the Π * arrangement in which a core of invariant tori is located on ( ∕2, 0, z) , a 2 1 screw axis.
The simplest field of all (206.4) klm,S in Table 3 is It has no remarkable invariant torus except very thin crooked tori. Indeed, Galloway and Proctor [29] reported that almost all integral curves are chaotic in numerically solving in which interchanging y and z on the left and right sides yields (62) As far as we know, the best example of the Π * arrangement of invariant tori with straight cores is observed in = (206.4) 001,S + (206.4) 021,S . Here, is obtained from Table 3. Because of the necessary condition (7) for the line ( ∕2, 0, z) to be a core of invariant tori reads that is, Figure 14 shows Π * -arranged invariant tori in the sum of (62) and (63) with Although obscure in the figure, the invariant torus along the line ( ∕2, 0, z) (resp. (3 ∕2, 0, z) ) is left-handed helically (resp. right-handed helically) wound around by integral curves. The volume percentage of the Π * -arranged ordered region is d dt   with (62), (63) and (64). One of the initial points is (10) with (n x , n y , n z ) = (10, 1, 1) estimated at 5.8% . This is not to be compared with the volume percentage of rods in the Π * rod packing ( 3 ∕16 = 58.90 … % , see [7,9]).
The Beltrami flow (214.3) 013,RB also has chain-like invariant tori, as shown in Fig. 15b. Each of the invariant tori, which is unknotted (although looking as if knotted in the figure), interweaves with another and forms a (12, 2)-torus link. In addition, this pair of invariant tori links with three other pairs. The volume percentage of the chain is estimated at 0.3%.
In [21, Fig. 8], we showed the − Ω arrangement of invariant tori with curved cores in The cores of the invariant tori are straightened by adding a 1 5 -wavelength field 013,RB with (71)-(73). a The + Ω arrangement when 1 and 2 are from (74). One of the initial points is (10) with (n x , n y , n z ) = (17, 29, 1) . b The − Ω arrangement when 1 and 2 are from (75). One of the initial points is (10) with (n x , n y , n z ) = (35, 18, 1) . c The I-WPlike arrangement when 1 and 2 are from (74). It is covered with a single integral curve that starts from (10) with (n x , n y , n z ) = (19,20,14) and are plotted modulo 2 Figure 17a and b depict the + Ω and − Ω arrangements of invariant tori with straight cores in the linear combinations with (74) and (75), respectively. The volume percentages of the + Ω -and − Ω-arranged ordered regions are estimated at 12.0% and 29.3% , respectively. By contrast, the percentage of rods in the + Ω or − Ω rod packing is √ 3 ∕18 = 30.22 … % [7]. In addition to the + Ω -and − Ω-arranged invariant tori, other invariant tori, whose genus is not equal to 1, are observed in the linear combination with (74) or (75). As shown in Fig. 17c, they are surfaces resembling Schoen's I-WP minimal surface, whose genus is 4 in a primitive kaleidoscopic tetrahedron [31] and 7 in a minimal cubic unit cell [32]. Unlike a usual invariant torus of genus 1, each of the I-WP-like invariant tori contains zeros of the vector field on it. The zeros are located at the points where two of x, y and z coordinates are i 1 and i 2 with integers i 1 and i 2 such that i 1 ≡ i 2 mod 2 . The volume percentage of the region occupied by the I-WP-like invariant tori has not been measured because some integral curves on them are difficult to distinguish from chaotic integral curves.

Summary and discussion
The sets of invariant tori whose arrangements resemble the cubic rod packings in Figs. 1-6 have been considered in solenoidal vector fields. In principle, they are obtained by the process in Sect. 2.2. For chiral arrangements, we used Beltrami flows (see (4)), which have the strongest chirality of all solenoidal vector fields satisfying (2). When the invariant tori are assumed to have straight cores, the sufficient condition (8) or (9) for their nonexistence is convenient to study which magnetic group to select in the process for each rod-packing arrangement. This answer is shown in Table 5. It was analytically derived from (8) or (9) with some symmetry or antisymmetry operations of associated magnetic groups or with the vector fields in Tables 1-4. If familiar with black/red-coloured symmetry/antisymmetry diagrams of magnetic groups in [26], we notice that the rotation or screw rotation axes that coincide with the cores of invariant tori are coloured black and orthogonally intersect with red-coloured twofold rotation axes. 5 The word "Maximal" in Table 5 implies that invariant tori may persist even if small perturbation that violates the symmetry or antisymmetry of the magnetic group is added to the vector field. A demonstration of Table 5 was given using solenoidal vector fields whose magnetic groups are listed in the table, as depicted in Figs. 10-17. If the cores of invariant tori are not assumed to be straight but allowed to curve, then another magnetic group may be selectable = 0.0446 … . 5 There are some colour misprints in twofold rotation axes in the diagrams of I P 4 1 32 and I P 4 ′ 1 32 ′ in [26].
for each arrangement. Even when restricted to vector fields of the magnetic groups in Table 5, a variety of invariant tori with curved cores are observed in addition to those with straight cores. Especially, chain-like arrangements of closed invariant tori (with or without knots) in Figs. 10b, c, 15b and 16b are notable. They correspond to the polycatenane lcv-y (see [5, Fig. 41]) or its variants. Furthermore, an ordered region formed by I-WP-like invariant tori, whose genus is not 1, is also observed as in Fig. 17c. In each vector field used for drawing Figs. 10-17, the volume percentage of the ordered region formed by nested invariant tori was estimated by counting the figures with invariant tori with the aid of supervised machine learning (based on the logistic regression method). There exist two kinds of errors in this measurement. One comes from the number 4N 3 of initial points (10). The larger it becomes, the more accurate the percentage becomes. The other comes from the error in the numerical solution (t) of (6). If the upper bound of integration t max is larger than our setting t max = 125 , then the error becomes larger and the ratio of figures with chaotic curves (resp. with invariant tori) goes up (resp. down). Judging from these factors, the volume percentages are considered to be reliable up to order 1%, although reported up to order 0.1%.
As mentioned in Sect. 1, some vector fields considered in this study can be realised in electromagnetic standing waves, in which invariant tori are wound around by electric or magnetic force lines (see [33, Fig. 1]). The application of these invariant tori to materials science may be an interesting topic.