Hedgehog-lattice spin texture in classical Heisenberg antiferromagnets on the breathing pyrochlore lattice

The hedgehog lattice, a three-dimensional periodic array of magnetic monopoles and antimonopoles, is known to be realized in the presence of the Dzyaloshinskii-Moriya (DM) interaction. Here, we demonstrate by means of Monte Carlo simulations that the hedgehog lattice is induced by not the DM interaction but frustration in classical Heisenberg antiferromagnets on the breathing pyrochlore lattice. In the model, the breathing bond-alternation is characterized by the ratio of the nearest-neighbor (NN) antiferromagnetic exchange interaction for large tetrahedra to that for small ones, J1'/J1. A quadruple-q state with the ordering vector of q=(\pm 1/2,\pm 1/2,\pm 1/2), which is realized for a large third NN antiferromagnetic interaction along the bond direction J3, turns out to become the hedgehog-lattice state in the breathing case of J1'/J1<1, while in the uniform case of J1'/J1 =1, it is a collinear state favored by thermal fluctuations. It is also found that in a magnetic field, the structure of the (1/2,1/2,1/2) hedgehog lattice is changed from cubic to tetragonal, resulting in a nonzero net spin chirality which in a metalic system, should yield a characteristic topological Hall effect.

The magnetic skyrmion and hedgehog are spin textures both characterized by the integer topological charge which corresponds to, in units of 4π, the total solid angle subtended by all the spins involved. The distinction between the two is in the real-space structure. In contrast to the skyrmion in two dimensions, the hedgehog extends over the three-dimensional real space, so that it has a sigular point in the texture center. In this respect, the hedgehog is sometimes called the magnetic monopole. Noting that the solid angle Ω ijk for three spins S i , S j , and S k is related with the scalar spin chirality χ ijk = S i · (S j × S k ) via Ω ijk = 2 tan −1 χ ijk 1+Si·Sj +Si·S k +Sj ·S k [46], the hedgehog can be understood as a topologically protected chirality object.
The hedgehog lattice is an alternating periodic array of the hedgehogs (monopoles) and anti-hedgehogs (antimonopoles) each having positive or negative nonzero scalar spin chirality. In other words, it is a long-range order (LRO) of the chirality. When viewed as a LRO of spin itself, the hedgehog lattice is a multiple-q state described by more than one ordering-wave-vectors q's [22,23,[25][26][27][28][30][31][32][33][34][35]. A prominent aspect of the hedgehog lattice is the Hall effect of spin chirality origin [47], i.e., the so-called topological Hall effect. In the non-centrosymmetric magnets having the DM interaction MnSi 1−x Ge x with 0.25 < x, the occurrence of the hedgehog lattice has been evidenced by the topological Hall effect together with the observation of multiple-q Bragg reflections [22][23][24][25][26][27][28]. Recently, the hedgehog lattice phase was found also in the DM-free centrosymmetric magnet SrFeO 3 [48,49], but the ordering mechanism is still unclear and no other candidate compounds have been reported so far. In this work, towards the exploration of new classes of magnets hosting the hedgehoglattice spin texture, we theoretically investigate effects of the frustration which in two dimensions, induces the skyrmion lattice [36]. A typical example of frustrated three-dimensional systems would be classical Heisenberg antiferromagnets on the pyrochlore lattice.
The pyrochlore lattice is a three-dimensional network of corner-sharing tetrahedra. In the presence of not only the nearest-neighbor(NN) antiferromagnetic exchange interaction J 1 but also the third NN one along the bond directions J 3 , it turns out that classical Heisenberg spins are ordered into a quadruple-q state with the four ordering vectors q = ( 1 2 , 1 2 , 1 2 ), (− 1 2 , 1 2 , 1 2 ), ( 1 2 , − 1 2 , 1 2 ), and ( 1 2 , 1 2 , − 1 2 ) in units of 2π a in the cubic basis [50][51][52] and that the spins are collinearly aligned due to the effect of thermal fluctuations. Since the above J 1 -J 3 pyrochlore antiferromagnet can potentially host the hedgehog lattice in the sense that the ordered phase is the multiple-q state, we try to suppress the spin collinearity. For this purpose, we consider the so-called breathing pyrochlore lattice consisting of an alternating array of small and large tetrahedra (see Fig. 1) [51,53].
In the breathing pyrochlore lattice, the bond alternation should be reflected in different values of the NN interactions on small and large tetrahedra, J 1 and J ′ 1 , so that its strength could be measured by the ratio J ′ 1 /J 1 . In the different context of the spin-lattice coupling, the effect of the breathing lattice structure has already been examined: within the collinear-spin manifold, the fully ordered quadruple-q state with q = (± 1 2 , ± 1 2 , ± 1 2 ) is changed into a partially disordered one for smaller values of J ′ 1 /J 1 [51]. This indicates that for the present Heisenberg spins, the thermally-driven collinear ( 1 2 , 1 2 , 1 2 ) state would get unstable with decreasing J ′ 1 /J 1 , possibly leading to the hedgehog lattice. As demonstrated in Fig.  1, this is actually the case. One can see from Fig. 1 that the quadruple-q ( 1 2 , 1 2 , 1 2 ) state on the breathing pyrochlore lattice is nothing but the hedgehog lattice (see Secs. II and III for details). This ( 1 2 , 1 2 , 1 2 ) hedgehog lattice appears as a result of the frustration in the absence of the DM interaction. In this paper, we further clarify that the hedgehog lattice in a magnetic field possesses a nonzero net spin chirality which should yield emergent phenomena such as the topological Hall effect.
This paper is organized as follows: In Sec. II, we introduce the spin Hamiltonian and discuss its ordering properties based on analytical calculations. The results of our MC simulations will be shown in Sec. III. In Sec. IV, we will discuss the effect of a magnetic field and demonstrate the emergence of the net spin chirality which results in an unconventional topological Hall effect quite different from the DM-induced one. The origin of this field-induced topological Hall effect will be addressed in Sec. V. We end the paper with summary and discussions in Sec. VI.

A. Model
The spin Hamiltonian we consider is given by where S i is the classical Heisenberg spin at the site i, S(L) , i, j , and i, j denote the summations over site pairs on the small (large) tetrahedra, the second NN pairs, and the third NN pairs along the bond direction, respectively. As explained in Sec. I, the breathing bondalternation is characterized by the ratio of the NN interaction for large tetrahedra to that for small ones, J ′ 1 /J 1 , and the third NN antiferromagnetic interaction along the FIG. 1: Hedgehog-lattice spin texture at zero field on the breathing pyrochlore lattice, where the spin configuration is obtained at T /J1 = 0.01 in the Monte Carlo (MC) simulation for the Hamiltonian (1) with J ′ 1 /J1 = 0.6, J2/J1 = 0 and J3/J1 = 0.7. Red, green, blue, and yellow arrows represent spins on the corners 1, 2, 3, and 4 of the small tetrahedra shown in the inset, respectively. The small tetrahedron at the center of the cubic unit cell is outlined by black dots. The upper and lower central small tetrahedra correspond to the hedgehog (monopole) and the anti-hedgehog (anti-monopole), respectively (see the main text in Sec. III). bond direction J 3 is crucial for the occurrence of the hedgehog lattice. Although the second NN interaction J 2 is not essential, we have included J 2 for completeness. When we discuss the in-field properties of the system in Sec. IV, we include the additional Zeeman term −H i S i,z in the Hamiltonian (1), where H is the intensity of a magnetic field. Throughout this paper, we take the cubic unit cell of edge length a = 1 shown in Fig. 1, and thus, the total number of spin N and the linear system size L is related via N = 16L 3 .
B. Effect of the breathing lattice structure on the spin ordering To get insight into roles of the breathing bondalternation, we first discuss ordering properties of this model in the mean-field (MF) approximation [54,55]. Introducing the Fourier transform S i = q S α q exp(iq · r i ) with the site index i = (α, r i ), we can rewrite H into the following form where the sublattice indices α = 1, 2, 3, 4 correspond to the four corners of a tetrahedron [see the insets of Figs. 1 and 2 (d)]. By analyzing the lowest eigen value of the 4-by-4 matrix J αβ (q), λ q , we find that for various sets of parameters we searched, possible ordering wavevectors are 2π a (± 1 2 , ± 1 2 , ± 1 2 ), 2π a (±1, ±1, ±1), and incommensurate ones (see Fig. 7 The associated eigen vectors are given by with the normalization factor c N = 1/ √ 1 + 3ε 2 and the dimensionless coefficient ε is zero for the uniform pyrochlore lattice with J ′ 1 /J 1 = 1, and gradually increases with increasing the breathing alternation (decreasing J ′ 1 /J 1 ). Note that J 2 does not show up in λ Qµ and U Qµ , and thus, irrelevant for the ordering properties.
In the ( 1 2 , 1 2 , 1 2 ) ordered phase, the mean field S α q is described as where U α Qµ represents the αth component of U Qµ . In the uniform case of J ′ 1 /J 1 = 1 and thereby ε = 0, all the U Qµ 's are orthogonal to one another, so that Φ Q1 , Φ Q2 , Φ Q3 , and Φ Q4 correspond to the polarization vectors of the spins belonging to sublattice 2, 3, 4, and 1, respectively. In other words, to make spins exist at each sublattice, the ordered phase should be the quadruple-q state involving all the Φ Qµ 's, and thus, other multiple-q states such as single-q, double-q, and triple-q states are not realized in the present system. Then, the Ginzburg-Landau (GL) free energy in the form expanded with respect to Φ Qµ is obtained as F GL /(N/4) = f 2 + f 4 + δf 4 with Only in the breathing case of J ′ 1 /J 1 < 1, ε is nonzero and thus, δf 4 is active. The most stable spin configuration is obtained by minimizing F GL under the constraint Noting that in the uniform case, each Φ Qµ determines the spin orientation axisP µ at each sublattice, we assume that even in the breathing case, Φ Qµ takes the form of Φ Qµ = S √ 8 e iθµP µ . Then, we have In the uniform case of J ′ 1 /J 1 = 1, δf 4 = 0, so that the relation among different θ µ andP µ cannot be determined. In other words, spins at different sublattices are uncorrelated at the MF level. As will be shown below, however, thermal fluctuations induce the collinear intersublattice correlation. In the breathing case of J ′ 1 /J 1 < 1, δf 4 is active and the minimization condition for δf 4 is θ 1 −θ 2 −θ 3 −θ 4 = π/2 andP µ ·P ν = −1/3 (µ = ν), suggestive of a nonzero value of the chirality χ ijk = S i ·(S j ×S k ) formed by any three spins on each tetrahedron. It will be shown in the following numerical calculations that this quadruple-q ( 1 2 , 1 2 , 1 2 ) state with the nonzero chirality is nothing but the hedgehog lattice. It should be emphasized that at least at the MF level, an infinitesimally small breathing alternation yields nonzero ε and resultantly, can induce the ( 1 2 , 1 2 , 1 2 ) hedgehog lattice.

III. RESULT OF THE MC SIMULATION AT ZERO FIELD
Now that the role of the breathing bond-alternation has been understood, we shall perform a rigorous analysis of the finite-temperature properties of the Hamiltonian (1) by means of Monte Carlo (MC) simulations. In this section, we consider the zero-field case of H = 0.
In our MC simulations, 2×10 5 sweeps are carried out under the periodic boundary condition and the first half is discarded for thermalization, where our 1 MC sweep consists of 1 heatbath sweep and successive 10 overrelaxation sweeps.
Observations are done at every MC sweep and the statistical average is taken over 4 independent runs starting from different random initial configurations. We identified the low-temperature ordered phase by observing various physical quantities such as the spin collinearity P = 3 2 1 N 2 i,j S i ·S j 2 − 1 3 and the spin and chirality struc- where the scalar chirality of the lth tetrahedron with its center-of-mass position R l is defined by χ(R l ) = i,j,k∈l th tetra χ ijk and the order of i, j, and k is defined in the anticlockwise direction with respect to the normal vector of the triangle formed by the three sites i, j, and k,n ijk , pointing outward from R l . In the same manner, we define the solid angle subtended by the four spins on the tetrahedron as Ω(R l ) = i,j,k∈l th tetra Ω ijk .
The temperature dependence of the specific heat C and the spin collinearity P obtained in the MC simulations for J 2 /J 1 = 0 and J 3 /J 1 = 0.7 is shown in Fig. 2 (a). In the uniform pyrochlore lattice with J ′ 1 /J 1 = 1, P starts developing at the transition temperature indicated by the sharp peak in C, and approaches the maximum value of 1 at the lowest temperature. Since the Hamiltonian (1) itself does not have an interaction to collinearize spins, thermal fluctuations should be relevant to the spin collinearity as is often the case with frustrated magnets [56,57]. With increasing the breathing bond-alternation, i.e., decreasing J ′ 1 /J 1 , such a temperature effect is gradually suppressed. At J ′ 1 /J 1 = 0.6, P develops on cooling similarly to the case of J ′ 1 /J 1 = 1, but it finally goes to zero at T = 0 via an intermediate phase explained below, suggesting the occurrence of a noncollinear ground state. At the smaller value of J ′ 1 /J 1 = 0.2, the ordered phase is completely masked by the noncollinear state with P = 0. To look into the details of the ordered phases, we examine the spin and chirality structure factors F S (q) and F C (q). Figure 2 (b) shows F S (q) in the noncollinear state with P = 0. This state is characterized by the magnetic Bragg peaks with the same intensity at the cubicsymmetric families of (± 1 2 , ± 1 2 , ± 1 2 ), so that it is a cubicsymmetric quadruple-q state. Such a quadruple-( 1 2 , 1 2 , 1 2 ) spin structure can commonly be seen irrespective of the value of J ′ 1 /J 1 . One can see from Fig 2 ) = 2 h,k,l=±1/2 F S (h, k, l) serves as the order parameter to characterize the transition from the paramagnetic phase in all the three cases of J ′ 1 /J 1 = 1, 0.6, and 0.2. Although whether the spins are collinear or noncollinear is indistinguishable in the spin sector F S (q), it is clearly visible in the chirality sector F C (q). In the uniform case of J ′ 1 /J 1 = 1, spins are collinearly ordered, so that the local scalar chirality χ ijk and the associated F C (q) are vanishingly small. In contrast, in the breathing case of J ′ 1 /J 1 < 1, F C (q) in the noncollinear state with P = 0 exhibits Bragg peaks at (± 1 2 , ± 1 2 , ± 1 2 ) suggestive of a tetrahedron-based chirality order [see for the noncollinear ( 1 2 , 1 2 , 1 2 ) state with P = 0. The real-space spin configuration of the noncollinear state is shown in Fig. 1. Spins belonging to each sublattice, which are represented by the same colored arrows in Fig. 1, constitute up-down antiferromagnetic chains along the bond directions, and four sublattices have different polarization axes. In Fig. 1, since four outgoing (incoming) spins on a tetrahedron subtend the solid angle of Ω(R l ) = 4π (−4π), the upper (lower) small tetrahedron at the center of the cubic unit cell has the topological charge of +1 (−1). As the spins are ordered in a way such that the two cubic unit cells are alternately arranged, the state is nothing but the hedgehog lattice in which the same number of the monopoles with +1 charge and the antimonopoles with −1 charge are equally populated and their densities n + and n − are given by n + = n − = 1 16 . For the monopole and antimonopole tetrahedra, the total spin on each tetrahedron is zero as is expected from the NN antiferromagnetic interaction J 1 or J ′ 1 [54], while not for other tetrahedra. This is because the ( 1 2 , 1 2 , 1 2 ) spin ordering in the present system is driven by the relatively strong third NN antiferromagnetic interaction J 3 rather than J 1 or J ′ 1 and thus, not all the tetrahedra satisfy the constraint of total spin zero stemming from J 1 and J ′ 1 . In relation to this, as one can see from Fig. 1, the monopoles and antimonopoles, which satisfy the constraint, are formed not on the large tetrahedra but on the small tetrahedra, since the latter has the stronger NN antiferromagnetic interaction J 1 > J ′ 1 and thus, tries to follow the constraint as much as possible. We note that in the present system, the isotropic Heisenberg spins spontaneously form this chirality order as a result of the frustration, which is in sharp contrast to spinice systems with strong local magnetic anisotropy where a similar ( 1 2 , 1 2 , 1 2 ) spin structure has been discussed [58]. In the quadruple-q ( 1 2 , 1 2 , 1 2 ) phases other than the hedgehog lattice, relative angles between the polarization axes for the four sublattices are changed with each updown antiferromagnetic chain being almost unchanged. Figure 2 (d) shows the map of the spins in the three different ( 1 2 , 1 2 , 1 2 ) phases obtained at J ′ 1 /J 1 = 0.6. One can see from the bottom (top) panel of Fig. 2 (d) that in the higher-temperature collinear (lower-temperature hedgehog-lattice) phase, the polarization axes for the four sublattices are the same (different). The spin ordering pattern in the collinear phase is up-up-down-down along all the bond directions [50,51]. As shown in the middle of Fig. 2 (d), the intermediate-temperature phase between the two is a coplanar state in which four sublattices are divided into two pairs and spins belonging to each pair have the same polarization axis (see Fig. 8 in Appendix B for details of the real-space spin configurations).
The stability region of the ( 1 2 , 1 2 , 1 2 ) hedgehog lattice is summarized in Fig. 3. Among the three quadrupleq ( 1 2 , 1 2 , 1 2 ) states, only the coplanar one cannot survive down to the lowest temperature, so that it does not appear in Fig. 3. The ( 1 2 , 1 2 , 1 2 ) hedgehog lattice is realized in the wide range of the parameter space even in the presence of J 2 . As expected from the MF analysis, a small breathing alternation, i.e., small deviation 1 − J ′ 1 /J 1 , seems to be sufficient to stabilize the hedgehog lattice.

IV. EFFECT OF A MAGNETIC FIELD
In this section, we will discuss the effect of a magnetic field H on the ( 1 2 , 1 2 , 1 2 ) hedgehog lattice. Figure 4 (a) shows the temperature and magnetic-field phase diagram for J ′ 1 /J 1 = 0.2, J 2 /J 1 = 0, and J 3 /J 1 = 0.7. The low-field and high-field phases are the ( 1 2 , 1 2 , 1 2 ) hedgehoglattice and coplanar states being subject to spin canting, respectively. We note that for larger values of J ′ 1 /J 1 , the canted collinear state also shows up (see Fig. 9 in Appendix C). Although a weak transition-like anomaly is found in the high-field phase [see open symbols in Fig.  4 (a)], the detailed analysis of the high-field phase is beyond the scope of this work, and hereafter, we will focus on the low-field hedgehog-lattice phase. As shown in the left panels of Fig. 4 (b), the in-field hedgehog lattice is characterized by the enhanced (suppressed) (± 1 2 , ± 1 2 , ± 1 2 ) magnetic Bragg reflections for the spin component parallel (perpendicular) to the field O Also, the magnetization m shows a weak downward convex curve as a function of H in this lowfield phase. Besides, in a magnetic field, additional weak Bragg peaks emerge at the wave vectors of the (0, 0, 1) and (1, 1, 0) families in both F S (q) and F C (q), which can clearly be seen in Fig. 4 (c). This in-field hedgehog lattice is tetragonal symmetric in the sense that among the three cubic symmetric points of (±1, 0, 0), (0, ±1, 0), and (0, 0, ±1), only one is selected, and such a situation is also the case for the (1, 1, 0) family. In the case of Fig.  4 (c), z-direction is special. In the chirality sector [the lower panels of Fig. 4 (c)], only (0, 0, 1) and (1, 1, 0) are picked up, whereas in the spin sector (the upper ones), the rest four, i.e., (1, 0, 0), (0, 1, 0), (1, 0, 1) and (0, 1, 1), are selected. As one can see from the right panels of Fig.  4 (b), the averaged intensities of these additional peaks in the spin and chirality sectors O S ⊥ (001) , O S ⊥ (110) , O C (001) , and O C (110) are nonzero only in the in-field hedgehog lattice.
When the system is metallic and localized spins S i are coupled to conduction electrons, the field-induced symmetry reduction from cubic to tetragonal in the hedgehog lattice is reflected in the Hall effect of spin chirality origin, i.e., the topological Hall effect. In the weakcoupling theory [47], the topological Hall conductivity σ T µν is proportional to the total spin chirality multiplied by a geometrical factor χ T = i,j,k S,L χ ijknijk , namely, σ T µν ∝ ǫ µνρ χ T ρ , where i, j, k S,L denotes the summation over all the three sites on each tetrahedron and n ijk represents the surface normal. χ T corresponds to the emergent fictitious magnetic field. Although according to Ref. [47], χ T should be calculated by taking further NN sites into account, we have restricted only to the NN triads because the dominant contribution comes from such short-distance triads. Indeed, we have checked that the original and present χ T 's exhibit qualitatively the same field dependence. One can see from the bottom right panel of Fig. 4 (b) that with increasing the external magnetic field, the fictitious field |χ T | is gradually induced in the hedgehog-lattice phase and it vanishes out of this topological phase. It is also found that the direction of χ T is not arbitrary but rather takes one of the six degenerate directions, ±x, ±ŷ, and ±ẑ. As inferred from the similar field dependences of 110) , and |χ T | [see right panels of Fig. 4 (b)], the preferred direction of χ T is determined from the tetragonal symmetry of the spin structure. For example, when the z-direction is picked up for the tetragonal spin structure like in the case of Fig. 4 (c), χ T is along theẑ-direction, so that σ T xy = −σ T yx takes a finite value while other components are zero. In the next section, we will address the origin of this unique type of the topological Hall effect characteristic of the ( 1 2 , 1 2 , 1 2 ) hedgehog lattice.

V. ORIGIN OF THE FIELD-INDUCED TOPOLOGICAL HALL EFFECT
To examine the association between the field-induced tetragonal-symmetric hedgehog-lattice and the net chirality χ T , or equivalently, the emergent fictitious field for the topological Hall effect, we shall first discuss the magnetic structure of the in-field hedgehog lattice. Figure 5 (a) shows the spin configuration associated with Fig. 4 (c). Since the ( 1 2 , 1 2 , 1 2 ) hedgehog lattice consists of the alternating array of two cubic unit cells having oppositely oriented spins [see Fig. 1], only one of the two is shown in Fig. 5. One can see from the lower panels of Fig. 5 (a) that the sublattices 1 and 2, and 3 and 4 are paired up and two spins in each pair are collinearly aligned in the S x S y spin plane. Due to this pairing, the symmetry of each tetrahedron is reduced to tetragonal with respect to  Once a pairing symmetry is fixed, there exist two types of the tetrahedral spin configurations, which we call AF-and F-types. As shown in Fig. 5 (a), in the AF-type (F-type) tatrahedron, the S x S y components in each pair are antiferromagnetically (ferromagnetically) aligned, whereas the S z ones ferromagnetically (antiferromagnetically). Also, the collinear axes for the two pairs are orthogonal to each other in the AF-type, while not in the F-type. The AF-type is further categorized into two: one involves S z 's of the same sign [tetra I in Fig. 5 (a)] and the other involves S z 's of opposite signs [tetra II in Fig. 5 (a)]. The monopole and antimonopole correspond to the latter. The distribution of these AF-and F-type tetrahedra within the cubic unit cell is also tetragonalsymmetric. One can see from Fig. 5 (d) that the AF-type tetrahedron pair and the F-type one are stacking along the tetragonal-symmetric z-axis. Thus, the symmetry-reduction to tetragonal occurs at the level of the cubic unit cell as well as each tetrahedron. Now that the spin structure of the in-field hedgehog lattice is clarified, we shall next discuss how it is related to the total spin chirality χ T = i,j,k S,L χ ijknijk . As a first step, we consider the associated local chirality χ T,tetra = i,j,k∈tetra χ ijknijk for the small tetrahedron shown in Fig. 6 (a), which is calculated as When the tetrahedron is projected onto the twodimensional planes perpendicular to the x-, y-, and zaxes [see Figs. 6 (b)-(d)], the chirality defined on each two-dimensional plane χ 2D,tetra is given by a summation of four different χ ijk 's with triad (i, j, k) ordered in the anticlockwise direction with respect to the surface normal. Then, one finds that the so-obtained χ 2D,tetra 's correspond to the x-, y-, and z-components of χ T,tetra in Eq. (8), which suggests that the ρ component of χ T is determined from the spin configurations on the projected two-dimensional planes perpendicular to the ρ-direction. Thus, we next consider how the tetragonal symmetry in each tetrahedron, i.e., the pairing pattern, looks on the projected planes.
One can see from wavy lines in Figs. 6 (a)-(d) that the paired spins are arranged diagonally on the plane perpendicular to the tetragonal-symmetric direction (zdirection in Fig. 6), while adjacently on other planes. Such a difference in the pairing arrangement yields an anisotropy in the local chirality. It turns out from a simple calculation (see Appendix D) that χ 2D,tetra on the projected plane becomes nonzero for the diagonal arrangement, while it is basically zero for the adjacent ones except a particular case of the F-type spin configuration, which corresponds to Fig. 10 Bearing the association between the tetragonal symmetry and the local chirality in our mind, we will turn to the total chirality on the projected two-dimensional planes, i.e. layers of one-tetrahedron width [see Figs. 6 (e) and (f)].
Figures 6 (g) and (h) show the layer-resolved distributions of χ T,tetra calculated from the spin configuration in Fig. 5. One can see from Fig. 6 (g) that on the layers perpendicular to the non-tetragonal-symmetric x-direction where the paired spins are arranged adjacently, the local chirality χ T,tetra x for most of the tetrahedra takes a vanishingly small value, which may be positive or negative depending on the thermal noise. Although an exceptional F-type tetrahedron, e.g., tetra III in the zoomed view in Fig. 6 (g), yields a nonzero large χ T,tetra x , it is completely canceled out by the contribution from the counter tetrahedron in the neighboring cubic unit cell, so that the net chirality vanishes on any yz-plane layer. Such a situation is also the case for the layers perpendicular to the y-direction. By contrast, on the layers perpendicular to the tetragonal-symmetric z-direction where the paired spins are arranged diagonally, not only the local chirality χ T,tetra z but also the net chirality becomes nonzero. As shown in Fig. 6 (h), χ T,tetra z 's on the layer without monopoles [layer 1 xy in Fig. 6 (h)] take nonzero values of the same sign, leading to a large total chirality, while χ T,tetra z 's on the layer with monopoles [layer 2 xy in Fig. 6 (h)] take positive and negative signs. As one can see from the zoomed view focusing on the two cubic unit cells in Fig. 6 (h), the monopole and antimonopole tetrahedra yield χ T,tetra z 's of the same sign and their contributions are canceled out by those from other tetrahedra. Thus, a dominant contribution comes from the layers without the monopoles, which has actually been confirmed from the numerical data for Fig. 6 (h) together with the result of analytical calculations shown in Appendix D.

VI. SUMMARY AND DISCUSSION
In this paper, we have demonstrated by means of MC simulations that the hedgehog lattice characterized by the ( 1 2 , 1 2 , 1 2 ) magnetic Bragg reflections is induced by the frustration in classical Heisenberg antiferromagnets on the breathing pyrochlore lattice, where the breathing bond-alternation quantified by the ratio of the NN antiferromagnetic interactions J ′ 1 /J 1 and a large third NN one along the bond directions J 3 are essential. It is also found that in a magnetic field, the structure of the ( 1 2 , 1 2 , 1 2 ) hedgehog lattice is changed from cubic to tetragonal, resulting in the nonzero net spin chirality χ T , or equivalently, nonzero emergent fictitious magnetic field, along the tetragonal-symmetric direction.
In the DM-induced hedgehog-lattice phase, the fictitious field is also caused by the magnetic field, but its origin is the field-induced position shifts of monopoles and antimonopoles [33], which is in sharp contrast to the present frustrated system where their positions are unchanged. As a result, the total spin chirality χ T appears in any of the possible three directions x, y, and z in the present system, while only along the applied field direction in the DM system. In addition, the sign of χ T can be positive or negative in the former, while it is fixed by the DM interaction in the latter. Such a degeneracy of the right-handed and left-handed chiralities is a unique aspect of the frustration-induced topological spin textures distinct from the DM-induced ones [36].
In general, an ordering wave-vector minimizing the eigen value of J αβ (q) defined in Eq. (2) should be realized in the ground state, provided that the local spinlength constraints |S i | = 1 are satisfied. Thus, searching the ordering wave-vector with the lowest eigen value is a starting point to analyze ordering properties of the system. Figure 7 shows ordering wave-vectors with the lowest eigen value of J αβ (q). For larger values of J 3 /J 1 , the ground-state candidate is a state with the ordering wave-vector of (± 1 2 , ± 1 2 , ± 1 2 ), whereas for smaller values of J 3 /J 1 , it is a state with an incommensurate wavevector. With increasing J 2 /J 1 or decreasing J ′ 1 /J 1 , the ( 1 2 , 1 2 , 1 2 ) state gets more stable. A commensurate (1, 1, 1) state is also possible for an antiferromagnetic J 2 and a weak J 3 .
Appendix B: Real-space spin configurations of the three different quadruple-q ( 1 2 , 1 2 , 1 2 ) states In the MC simulations, we obtain the three different quadruple-q ( 1 2 , 1 2 , 1 2 ) states, the collinear, coplanar, and hedgehog-lattice ones. Their real-space spin configurations at H = 0 are shown in Fig. 8. In all the three cases, spins belonging to each sublattice (the same colored arrows in Fig. 8) constitute up-down antiferromagnetic chains along the bond directions, although exactly speaking, the spin directions in each chain are not perfectly collinear when the lattice is not uniform, which will be explained below.
The collinear phase can survive down to a low temperature basically for the uniform pyrochlore lattice with J ′ 1 /J 1 = 1. Since the polarization axes of the different sublattices are collinearly aligned, the resultant spin configuration consists of the up-up-down-down antiferromagnetic chains running along all the bond directions. In this case, the spin directions in each chain are perfectly collinear.
The coplanar phase is realized only at moderate temperatures. Indeed, when we evaluate the free energy of the coplanar state by using Eq. (7), it is higher than that of the hedgehog lattice. If the coefficient A 3 defined in Eq. (6) were large enough, the coplanar state could be stable, but this is not the case in the present model. Thus, higher-order contributions in the GL free energy, which have not been considered in Eq. (7), or thermal fluctuations may be relevant for the stability of the coplanar phase. One can see from Fig. 8 (b) that in the coplanar phase, four sublattices represented by different colors are divided into two pairs and spins belonging to each pair have the same polarization axis, suggesting that the spin state is coplanar although the coplanarity is disturbed by the thermal noise. We note that as we will discuss below, even if this state were realized at a low temperature escaping from the thermal noise, the spins in each pair could not be perfectly collinear.
The hedgehog-lattice phase is realized at low temperatures in the wide range of the parameter space. Although the real-space spin configuration at J ′ 1 /J 1 = 0.6 and T /J 1 = 0.01 is presented in Fig. 1, here, we show the one for the stronger breathing alternation of J ′ 1 /J 1 = 0.2 in Fig. 8 (c). The snapshot is obtained by cooling the system down to the very low temperature T /J 1 = 0.001 to reduce the thermal noise. One can see from the lower panel of Fig. 8 (c) that the polarization vector of each sublattice splits into four, although the real-space spin configuration in the upper panel does not differ so much from Fig. 1. As we will address below, such a splitting is generic to the quadruple-q ( 1 2 , 1 2 , 1 2 ) noncollinear states on the breathing pyrochlore lattice. Actually, when the spin configuration in Fig. 1 is further cooled down, the splitting overcomes the thermal noise and becomes visible. The origin of the splitting can be understood within the MF approximation.
Since the Fourier component of spin S α q is obtained with the combined use of Φ Qµ = S √ 8 e iθµP µ and Eqs. (3) and (5), spins belonging to sublattice α, S α j 's, are given by +P 3 cos(Q 3 · r j + θ 3 ) +P 4 cos(Q 4 · r j + θ 4 ) , As long asP µ 's are not collinear like in the cases of the coplanar and hedgehog-lattice phases, spins on each sublattice cannot be collinear because of the existence of ε. Particularly in the hedgehog-lattice phase, all the four spins on each sublattice within the cubic unit cell orient in the different directions, exhibiting the four spots when the spins are mapped onto the sphere. Since ε is nonzero only in the breathing case, such a splitting is inherent to the ( 1 2 , 1 2 , 1 2 ) noncollinear states on the breathing pyrochlore lattice. We note that ε is important for the occurrence of the hedgehog-lattice phase as explained in Sec. II, but the resultant splitting itself is not. Appendix C: The temperature and magnetic-field phase diagram in the weaker breathing case of J ′ 1 /J1 = 0.6 Figure 9 shows the temperature and magnetic-field phase diagram for J ′ 1 /J 1 = 0.6, J 2 /J 1 = 0, and J 3 /J 1 = 0.7. Compared with the strongly breathing Fig. 4 (a)], the canted collinear state additionally shows up in the higher temperature region, and the stability region of the hedgehog lattice is suppressed. In the canted coplanar state, a weak transition-like anomaly is also found for J ′ 1 /J 1 = 0.6 as well as J ′ 1 /J 1 = 0.2 (see open symbols in Fig. 9).
Appendix D: Analytical calculation of the local and total chiralities in the in-field hedgehog-lattice phase As shown in Figs. 6 (a)-(d), when a tetragonalsymmetric tetrahedron is projected onto a twodimensional plane, two spins in each pair are arranged adjacently or diagonally on the plane. Here, we analytically calculate the chirality on the two-dimensional plane for various spin arrangements of a single tetrahedron shown in Fig. 10, and apply the result on the single tetrahedron to the in-field hedgehog lattice involving the 8 small tetrahedra.
We first consider the AF-type tetrahedron in which the S x S y components in each pair are antiferromagnetically aligned and the S z ones ferromagnetically. Since the scalar chirality is invariant under any rotations with respect to the applied field direction, i.e., S z , we can fix a S x S y component of one spin to be, for example, in the S y direction. Thus, we assume that two spins in one pair are given by Then, the rest two spins in the other pair can be expressed as where φ determines the angle between the S x S y -plane collinear axes of the two pairs. The single-tetrahedron chirality on the two-dimensional plane is given by χ 2D,tetra = i,j,k χ ijk , where i, j, k denotes the summation over all the four triads and χ ijk = S i · (S j × S k ) is defined in a way such that i, j, and k are ordered in the anticlockwise direction with respect to the surface normal. For the spin configurations shown in Figs. 10 (a)-(d), the single-tetrahedron chirality χ 2D,tetra can be calculated as Fig. 9(a)] 0 [ Fig. 9 Fig. 9(d)].
(D4) Next, we consider the F-type tetrahedron in which the S x S y components in each pair are ferromagnetically aligned and the S z ones antiferromagnetically. and other notations are the same as those in Fig. 11. In each tetrahedron, the SxSy components of the two different pairs are assumed to be orthogonal to each other for simplicity.
Then, the rest two spins in the other pair can be expressed as S T2 = (S + xy cos φ, S + xy sin φ, m + δS z ) S T2 ′ = (S − xy cos φ, S − xy sin φ, m − δS z ). (D6) The single-tetrahedron chirality on the two-dimensional plane χ 2D,tetra can be calculated as (D7) As the chirality for a single tetrahedron χ 2D,tetra is obtained, we will turn to the total chirality in the hedgehog lattice phase. The spin configuration associated with Figs. 5 and 6 is shown in Fig. 11. Note that the spin structure is tetragonal-symmetric with respect to the zaxis. This spin configuration is schematically drawn in the layer-resolved form in Fig. 12, where for simplicity, φ is assumed to be 0 or π such that the angle between the S x S y -plane collinear axes of the two pairs be orthogonal to each other. Since the hedgehog lattice consists of the alternating array of the two cubic unit cells having oppositely oriented spins except the uniform magnetization m, the two neighboring unit cells on the projected plane (two neighboring unit cells in Fig. 12) involve the contribution given by Eqs. (D4) and (D7) and the counter contribution obtained by carrying out the replacement S ± xy → − S ∓ xy , δS z → − δS z (D8) in Eqs. (D4) and (D7). On the plane perpendicular to the non-tetragonalsymmetric direction [layers 1 yz and 2 yz in Fig. 12 (a)], the paired spins are arranged adjacently, so that χ 2D,tetra becomes zero for almost all the tetrahedra except tetra III, which corresponds to Fig. 10 (f), and its counter tetrahedron in the neighboring unit cell. The tetra III yields the nonzero χ 2D,tetra given by Eq. (D7) for Fig.  10 (f), but it is completely canceled out by the contribution from the counter tetrahedron because of Eq. (D8). Hence, the total chirality on the two-dimensional plane vanishes when the two paired spins are arranged adjacently.
In contrast, on the plane perpendicular to the tetragonal-symmetric direction [layers 1 xy and 2 xy in Fig.  12 (b)] where the two paired spins are arranged diagonally, the local chirality χ 2D,tetra becomes nonzero for all the tetrahedra. Since the four tetrahedra in the cubic unit cell [tetra's I, II, III, and IV in Fig. 12 (b)] are not equivalent to one another, we will take the δS z components for these tetrahedra, δS z,I , δS z,II , δS z,III , and δS z,IV , as independent variables. After averaging over all the tetrahedra within the two unit cells, we obtain the total chirality per tetrahedron as for the layer with the monopoles [layer 2 xy in Fig. 12 (b)], where S ± xy,λ (λ =I, II, III, and IV) is defined by S ± xy,λ = 1 − (m ± δS z,λ ) 2 . One can see from Eq. (D9) that the total chirality for the layer without the monopoles definitely takes a nonzero value. In Eq. (D9), the sign of the total chirality is positive since it is calculated in the specific case of Fig. 12 (b), but, in general, it can be positive or negative depending on the spin configuration. Actually, when the spin-space inversion S x → −S x is performed for all the spins in Fig. 12 (b), the state is turned to the other chirality-degenerate state with its energy unchanged because the spin Hamiltonian is invariant under the spin-space inversion, and then, the total chirality changes its sign. On the other hand, on the layer with the monopoles [see Eq. (D10)], the cancellation occurs between the different kinds of tetrahedra, and as a result, the net chirality on this layer becomes small. Thus, the overall chirality is dominated by the contribution from the layers without the monopoles, which has actually been confirmed from the numerical data calculated from the MC spin snapshots. Since as one can see from Eqs. (D9) and (D10), the total chirality becomes non-vanishing only for nonzero m, the magnetic field as well as the symmetry reduction to tetragonal is crucial for the emergence of the topological Hall effect in the present system.