{'fraction_non_alphanumeric': 0.04972828995089004, 'fraction_numerical': 0.01809814856653626, 'mean_word_length': 4.003450737743932, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 28, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Dynamical Piezomagnetic Effect in Time-Reversal Invariant Weyl Semimetals with Axionic Charge-Density Waves

Charge-density waves (CDWs) in Weyl semimetals (WSMs) have been shown to induce an exotic axionic insulating phase in which the sliding mode (phason) of the CDW acts as a dynamical axion field, giving rise to a large positive magneto-conductance [Wang et al., Phys. Rev. B 87, 161107(R) (2013); Roy et al., Phys. Rev. B 92, 125141 (2015); J. Gooth et al., Nature 575, 315 (2019)]. In this work, we predict that dynamical strain can induce a bulk orbital magnetization in time-reversal- (TR-) invariant WSMs that are gapped by a CDW. We term this effect the"dynamical piezomagnetic effect"(DPME). Unlike in [J. Gooth et al., Nature 575, 315 (2019)], the DPME introduced in this work occurs in a bulk-constant (i.e., static and spatially homogeneous in the bulk) CDW, and does not rely on fluctuations, such as a phason. By studying the low-energy effective theory and a minimal tight-binding (TB) model, we find that the DPME originates from an effective valley axion field that couples the electromagnetic gauge field with a strain-induced pseudo-gauge field. In particular, whereas the piezoelectric effects studied in previous works are characterized by 2D Berry curvature, the DPME represents the first example of a fundamentally 3D strain effect originating from the Chern-Simons 3-form. We further find that the DPME has a discontinuous change at a critical value of the phase of the CDW order parameter. We demonstrate that, when there is a jump in the DPME, the surface of the system undergoes a topological quantum phase transition (TQPT), while the bulk remains gapped. Hence, the DPME provides a bulk signature of the boundary TQPT in a TR-invariant Weyl-CDW.

The last four decades have witnessed a paradigm shift in condensed matter physics driven by the discovery of the geometric phase and topology of electronic wave functions [1,2]. The search for experimentally observable response signatures of bulk nontrivial topology has emerged as central to the advancement of solid-state physics and material science. In two-dimensional (2D) systems, the Berry curvature in the momentum space not only provides an essential contribution towards quantized topological response effects but also provides nonnegligible contributions to various non-quantized physical phenomena. As an example of a quantized topological response, integrating the Berry curvature over the 2D (magnetic) first Brillouin zone (1BZ) has been shown to yield a topological invariant called the first Chern number, which determines the quantized Hall conductance in the quantum Hall effect or the quantum anomalous Hall effect (QAHE) [3][4][5][6]. A non-quantized example is the intrinsic anomalous Hall effect in a variety of ferromagnetic metals, which is given by the integration of the Berry curvature over the occupied states [7]. Other non-quantized response effects include the expression for orbital magnetic moments in terms of the Berry curvature [8][9][10] and the nonlinear Hall effect induced from a Berry curvature dipole [11][12][13][14][15]. More recently, it has been demonstrated that the piezoelectric response can also be related to the Berry curvature [16][17][18][19][20], and can have a discontinuous change across a TQPT in 2D TR-invariant systems [21]. A natural question to ask is whether a strain-induced response can be related to the topological property of three-dimensional (3D) systems.
In 3D systems, a celebrated topological response originates from effective axion electrodynamics through the action [22,23] S ef f,θ = e 2 32π 2 where (and for the remainder of this work unless specified otherwise) we choose units in which = c = 1, and notation in which duplicated indices are summed over for notational simplicity. In Eq. (1), F µν = ∂ µ A ν − ∂ ν A µ is the field strength of the electromagnetic U (1) gauge field A µ , and θ is the effective axion field. The bulk average value of θ in a 3D gapped crystal with vanishing Hall conductivity, labeled as θ bulk , is determined by the Chern-Simons 3-form [22] instead of the 2D Berry curvature, and it is only well defined modulo 2π. The winding of θ bulk in a fourth dimension gives the second Chern number [22]. As θ bulk has the same transformation properties as E · B, any symmetry that flips the sign of E · B (i.e., "axion-odd" symmetries) can quantize θ bulk to 0 or π modulo 2π, providing a symmetry-protected Z 2 indicator of axionic bulk topology in 3D insulators [22][23][24][25][26][27][28][29][30][31][32][33][34][35].
Various axion-odd symmetries can also provide simplified expressions for evaluating θ bulk , which are consistent in a 3D insulator with vanishing Hall conductivity; thereby the mismatch among the expressions indicates the existence of the gapless points or a nonzero Hall conductivity [36]. A direct physical consequence of nonzero θ bulk is the magnetoelectric effect [22], where an external electric (magnetic) field induces a magnetization (polarization) in the parallel direction; if quantized by an axion-odd symmetry, the effect is called the topological magnetoelectric effect [22,37,38]. Other physical consequences of a non-zero θ bulk include the surface half QAHE [22], a giant magnetic-resonance-induced current [39][40][41], the topological magneto-optical effect [42][43][44] (especially its exact quantization [42]), the zero-Hall plateau state [45][46][47], and the image magnetic monopole [48].
In this work, we explore the role of a "valley-separated" variant of an effective axion field in the physical response induced by dynamical strain. In particular, we propose that dynamical homogeneous strain can induce a nonzero bulk-uniform magnetization in a class of 3D TR-invariant insulators that have vanishing total θ bulk (thus vanishing magnetoelectric response) but have a non-quantized and tunable θ per valley. We term the strain-induced magnetoelectric response the DPME. Specifically, for the DPME to be relevant to experimental probes, we require a 3D TR-invariant insulator to have (i) low-energy physics that are well captured by a pair of TR-related valleys in the 1BZ, and (ii) non-vanishing valley axion fields, despite exhibiting an overall trivial θ bulk . As we will demonstrate below, because the DPME originates from a valley axion field, the DPME is related to the bulk 3D topological properties, and in particular, a discontinuous change in the DPME serves as a direct signature of boundary TQPT.
Following Ref. [49][50][51], we recognize that the requirements for observing the DPME are satisfied by a TRinvariant WSM with a bulk-constant CDW. We emphasize that in this work, the CDW order parameter is taken to be constant (i.e. static and homogeneous) in the bulk of the system, in contrast to Ref. [49][50][51], which focused on bulk fluctuations of the CDW order parameter. By studying the low-energy effective theory of a TR-invariant Weyl-CDW, we find that the phase of the CDW order parameter determines the valley axion field, leading to effective valley axion electrodynamics, whereas the CDW wavevector induces a valley layered QAHE. Importantly, TR symmetry restricts that the valley axion field and valley layered QAHE sum to zero when taken over two TR-related valleys, as TR symmetry reverses the signs of θ and the Hall conductivity. On the other hand, dynamical strain can act as a pseudo-gauge field in WSMs [52][53][54][55][56][57][58][59][60][61][62][63], which was previously proposed as a means of studying the chiral anomaly of WSMs. In this work, we find that pseudo-gauge fields can further couple to the electromagnetic field and the (valley) axion field in the effective action in the presence of a CDW. According to the effective action, the analogues of the valley axion electrodynamics and the valley layered QAHE for the pseudo-gauge field are the DPME and the piezoelectric response, respectively, which are non-vanishing after summing over valleys.
Crucially, the bulk response coefficient of the DPME is the average value of valley axion field that is determined by the Chern-Simons 3-form. As the Chern-Simons 3form can only exist in three or higher dimensions, the DPME has a completely different origin compared to the previously studied piezoelectric effects, which instead are due to the 2D Berry curvature [16][17][18][19][20][21]. As a result, the DPME cannot be realized by trivially stacking of 2D systems with nontrivial piezoelectric effects. In this sense, the generalization from 2D piezoelectricity to the 3D DPME is analogous to that from the 2D quantum spin-Hall insulator to the 3D TR-invariant strong topological insulator. Specifically, the 3D TR-invariant strong topological insulator cannot be constructed from a simple 3D stacking of the 2D quantum spin-Hall insulators, and was instead discovered through more involved theoretical efforts [1,2]. Hence, the DPME cannot be viewed as a simple 3D generalization of 2D piezoelectricitythe DPME is instead a intrinsically 3D effect.
To test the low-energy results against the UV completion, we construct a minimal TB model for a TRinvariant WSM with varying values of the CDW order parameters. We specifically demonstrate that a CDW can drive a minimal TR-invariant WSM into a weak topological insulator (WTI) phase with a nontrivial weak Z 2 index, originating from the odd valley Chern number induced by the CDW wavevector. By evaluating the bulk average value of the valley axion field at various phases of the CDW order parameter, we numerically verify the relation between the phase of the CDW φ and the valley axion field θ predicted by the low-energy theory. We further show that dynamical strain generically induces a nonuniform current distribution, which can be split into a uniform background current and a nonuniform part. The uniform background current gives rise to the piezoelectric effect, which is a trivial stacking of the 2D piezoelectric effect discussed in Ref. [21]. Conversely, we find that the nonuniform part of the current is mainly localized at the surfaces, and possesses opposite signs on two surfaces perpendicular to the CDW wavevector (provided a large enough sample), which induces the bulk orbital magnetization that corresponds to the intrinsically 3D DPME.
In addition, we demonstrate the existence of discontinuous changes in the DPME by continuously varying the bulk CDW phase in the TB model. We show that the discontinuous change originates from gap closings on certain system boundaries (appearing as two TR-related 2D gapless Dirac cones on one surface), which reverses the direction of the corresponding surface currents and thus dramatically changes the bulk orbital magnetization. We note that, in this work, we employ the high-  energy convention in which 2D and 3D Dirac cones have two and four components, respectively. This convention is notably distinct from the recent condensed-matter works [64][65][66], in which both 2D and 3D fourfold degeneracies are termed "Dirac cones." We observe that the gap closing on each Weyl-CDW surface takes the same form as the 2D Z 2 TQPT (i.e., the TQPT between a 2D TR-protected Z 2 trivial insulator phase and a 2D TRprotected Z 2 nontrivial topological insulator (TI) phase) at generic momenta, and is unaccompanied by a bulk gap closing; hence the surface gap closing, in the presence of two valleys, is a boundary TQPT that changes the relative surface Z 2 index. When the gap closings on different surfaces happen simultaneously, we can attribute the DPME jump to a 2π change in the valley axion field induced by the surface Z 2 TQPT, through the bulk-boundary correspondence of the low-energy effective theory. Our findings suggest that jumps in the DPME can serve as bulk signatures of boundary TQPTs.
This work is organized as follows. First, we will present an intuitive picture for the DPME in Sec. II. Then, we will describe a low-energy theory for the DPME in Sec. III, which we will then verify with a TB model in Sec. IV. In Sec. V, we will elucidate the relationship between the jump of the DPME and the boundary TQPT. Finally, we will conclude in Sec. VI by introducing a proposal for measuring the DPME in experiment.

II. INTUITIVE PICTURE
Before presenting supporting analytic and numerical calculations, we will first provide an intuitive picture of the DPME, in comparison with the Berry-curvature contribution to the piezoelectric response [16][17][18][19][20][21]. We start with a 2D TR-invariant insulator whose low-energy physics is described by two 2D gapped Dirac cones (regarded as two valleys below), as shown in Fig. 1(a). Because the two valleys are related by TR symmetry, there must be an oppositely-signed valley Hall conductance within each valley (determined by the integral of the Berry curvature over each valley), and the Hall currents induced by a uniform electric field must point in opposite directions and exactly cancel. However, a nonvanishing charge current can still be generated by applying a pseudo-electric field that points in opposite direction in each valley, which causes the induced Hall currents to point in the same directions and add up to a nonzero total current ( Fig. 1(b)). The pseudo-electric field can be generated by a dynamical homogeneous strain tensor u as E pse ∼u [54,62], and the resultant total Hall current j ∼u characterizes the piezoelectric effect [16][17][18][19][20]. With regards to TR symmetry, a static external electric field preserves TR symmetry whereas a static external pseudo-electric field breaks TR symmetry, explaining the dramatic differences in the current response. Although the piezoelectric effect is not quantized in one insulating phase, its nonzero discontinuous change across a 2D TR-invariant TQPT (e.g. a 2D Z 2 TQPT [67]) is proportional to the change of a topological invariant [21] and thus serves as a new experimental probe of the topology (or more precisely a change in topology).
Although the strain also acts as a pseudo-electric field in the 3D DPME, the 3D DPME is fundamentally distinct from a trivial stacking of the 2D piezoelectric effect, as the 3D DPME originates from 3D axion field instead of the Hall conductance. To see this, consider a TRinvariant 3D gapped system whose low-energy physics is captured by two 3D gapped Dirac cones with complex masses (also regarded as two valleys below), as shown in Fig. 1(c) and (d). Because the TR symmetry maps one valley to the other, a single Dirac cone is not necessarily TR-invariant, and, in the absence of additional crystal symmetries, generically has an unpinned effective axion field θ. For a finite-sized configuration of the 3D system with a fully gapped TR-invariant 2D boundary, TR symmetry requires that the two Dirac cones have opposite θ angles, implying that an external electric field would induce opposite magnetizations in the two valleys ( Fig. 1(c)), which sum to zero. In contrast, a pseudoelectric field, which can be induced by dynamical homogeneous strain and points in opposite directions at each valley, would induce the same magnetizations at the two valleys ( Fig. 1(d)), resulting in a nonzero total magnetization M ∝ E pse ∼u. We term this effect the DPME. The DPME is different from the conventional piezomagnetic effect [68] because the magnetization of the for-mer is proportional to the time derivative of the strain tensoru, whereas the magnetization of the latter is directly proportional to the strain tensor u. Crucially, the bulk average value of the effective axion field θ is determined by the integrals of the Chern-Simons 3-form [22][23][24] where A is the non-Abelian Berry connection, and the Chern-Simons 3-form can only exist in three or higher dimensions. Therefore, unlike the 2D piezoelectric effects, the θ-induced 3D DPME originates from the Chern-Simons 3-form instead of the Berry curvature, meaning that the DPME is intrinsically 3D and cannot be given by trivial stacking systems with a 2D piezoelectric effect.

III. LOW-ENERGY EFFECTIVE THEORY OF TR-INVARIANT WSMS WITH AXIONIC CDWS
In this section, we will provide a low-energy theory for the DPME in TR-invariant WSMs with axionic CDWs. We would like to emphasize that the derivation below is not confined to Weyl-CDWs, and can be generalized to any TR-invariant system with valley axion fields.

A. Minimal Model
A WSM phase can only emerge in systems that either break TR symmetry (magnetic materials) or break inversion symmetry (non-centrosymmetric crystals). CDWs in magnetic WSMs have been previously studied in numerous works, including Ref. [49 and 50]. In this work, we focus on CDWs in TR-invariant WSMs, which can be realized in non-centrosymmetric crystals. Since two Weyl points related by TR symmetry share the same chirality, and because the total chirality of the whole system must vanish [69], then there must be four Weyl points in a minimal model of a TR-invariant WSM. Fig. 2(a) schematically shows a distribution of four Weyl points in a minimal TR-invariant WSM, where for simplicity, we have enforced an additional mirror symmetry that flips y, labeled as m y , such that all four Weyl points are symmetry-related. The momenta of the four Weyl points take the form where α = ± indicates the relative chirality of the Weyl points, and a = 1, 2 is termed the "valley index." TR symmetry, labeled as T , relates k 1,α to k 2,α with the same chirality index α, while the mirror m y changes the chirality of Weyl points and relates k 1,α to k 2,−α . We would like to emphasize that, although we keep m y in the following derivation for simplicity, mirror symmetry is not essential for the physics discussed below. Through a unitary transformation of the bases and by rotating and rescaling axes, we can always transform the low-energy Lagrangian of the four Weyl points into the following form [49] where ψ t,r,a,α is a two-component field for the two bands that form the Weyl point at k a,α , σ 0,x,y,z are the Pauli matrices, v i indicates the Fermi velocity along the i direction, and t and r are time and position, respectively. In this work, we will for simplicity focus on the case in which i = x, y, z are the three laboratory directions. Following the derivation in Ref. [49 and 70], we keep k = 0 as the momentum-space origin of all fermion fields, as this choice naturally includes the Weyl-point-induced valley Hall effect in the effective action, as discussed below. Throughout this section on the low-energy theory, we adopt a proper rescaling of the space and fields to cancel the Fermi velocities as elaborated in Appendix. C; the Fermi velocities will later be restored for comparison to the TB model in Sec. IV.
As mentioned above, both T and m y change the valley index a of the fields in Eq. (3), while m y (T ) changes (preserves) the chirality index α. Thus, we can always choose the bases to represent T and m y as iσ y K and −iσ y for the band index, respectively, where K is complex conjugation. According to the above symmetry representation for T and m y , a symmetry-preserving mean-field CDW term that couples two Weyl points of the same valley index a can be written as [49] L CDW = − a m a (r)e i(−1) a−1 Q·r ψ † a,+ ψ a,− + h.c. , (4) where the t, r dependence of ψ is implied, and where Q = k 1,+ − k 1,− = −(k 2,+ − k 2,− ) is the CDW wavevector, as shown in Fig. 2(a). Throughout this work, we will include the spatial dependence of the CDW order parameter m a (r), while keeping the order parameter timeindependent (i.e. static). TR symmetry requires that m 1 (r) = m 2 (r) * ≡ m(r), and m y symmetry requires that m * 1 (m y r) = m 2 (r). In general, m(r) = |m(r)|e iφ(r) is complex, and |m(r)| and φ(r) are the magnitude and phase of the CDW order parameter, respectively. As discussed in Sec. I, we consider the case where m(r) is equal to a complex constant mass m 0 = |m 0 |e iφ0 in the bulk throughout the work, i.e., |m(r)| = |m 0 | and φ(r) = φ 0 for r in the bulk. The underlying interaction that gives rise to the bulk CDW is discussed in Appendix. A at the mean-field level. Nevertheless, m(r), as well as |m(r)| and φ(r), can still have spatial dependence if the sample size is finite. Inspired by Ref. [22], we set |m(r)| → ∞ and φ(r) = 0 for r deep in the vacuum. Different gapped and symmetry-preserving boundaries can then be represented by different ways of smoothly connecting the bulk and vacuum limits of |m(r)| and φ(r). We next introduce the γ matrices where µ = 0, 1, 2, 3 and τ 0,x,y,z are Pauli matrices for the chirality index α. Using the above definitions of the γ matrices, we can rewrite the CDW term as where ψ a = ψ † a γ 0 , and Φ a (r) = (−1) a−1 (φ(r) + Q · r). For the remainder of this work, the spatial dependence of |m|, φ, Φ a will be implicit, and we will suppress the explicit dependencies on r for notational simplicity.
In order to elucidate the strain-induced linear response, we next introduce an electron-strain coupling for normal strain (i.e. stretch or compression along a specified axis) along the z direction, labeled as u zz (t). We require that the strain be adiabatic, homogeneous, and infinitesimal. Enforcing TR and mirror symmetries, the most general form of the leading-order electron-strain coupling reads where the time dependence of u zz is implied, and where the parameters ξ 0,x,y,z are material-dependent. In Eq. (7), we do not include the effects of strain that couple different Weyl points, as Weyl-point coupling strain is necessarily proportional to |m 0 |u zz , and because |m 0 | is typically small in real materials. The detailed procedure of adding the electron-strain coupling is shown in Appendix. B. We set the strain to be uniform throughout all of space, such that the gapped and symmetry-preserving boundary is implemented by the spatial dependence of the CDW order parameter m, as opposed to an inhomogeneous strain field.
(9) The valley-dependent U (1) gauge field takes the form which contains the physical gauge field A µ and the pseudo-gauge field induced by the strain u zz . In particular, the y component of the pseudo-gauge field can provide a pseudo-electric field that points in opposite directions in each of the two valleys E pse a = (−1) a ξ y eu zz e y .
As we will show below, all the nontrivial leading-order linear response comes from the pseudo-electric field in Eq. (11).

B. Effective Action
The low-energy response to A and u zz can be derived from the total effective action S ef f = a S ef f,a , which takes the form where the measure of the functional integral is in real space Dψ a Dψ a ∝ t,r dψ t,r,a dψ t,r,a .
We note that Eq. (3), as well as Eq. (12), are only exact when the momenta of the fermion fields in Eq. (13) are restricted near the Weyl points (or equivalently the momentum deviation from the corresponding Weyl point is below a finite momentum cutoff Λ). In Appendix. C, we demonstrate, however, that the correction to the response of interest brought by a finite Λ is of the order O(|m 0 | 2 /Λ 2 ), which is negligible owing to |m 0 | Λ in realistic materials. Since the focus of this work is on the low-energy response of the system, we can limit Λ → ∞ and take the functional integration (Eq. (13)) over the entire k ∈ R 3 , following Ref. [49-51, and 70].
The physical U (1) gauge field A and the strain tensor u zz in Eq. (8) are treated as fixed backgrounds, meaning that we neglect their dynamics. Under this assumption, L a has local valley U (1) gauge invariance, i.e., invariance under ψ a → ψ a e ieΓa(t,r) and A a,µ → A a,µ − ∂ µ Γ a (t, r), where Γ a is a valley-dependent scalar function and the corresponding transformation on ψ a is implicit here (and will remain implicit and for the reminder of this work). The valley U (1) gauge invariance corresponds to a separate vector current conservation for each valley; as different valleys are decoupled in Eq. (8), we preserve the current conservation in each valley against all orders of quantum correction, which is reasonable as long as valleys are well defined. As a result, the measure of the functional integral (Eq. (13)) is invariant under ψ a → e ieΓa ψ a , and S ef f,a is gauge in- µ ] allows us to perform a gauge transformation A a,µ → A a,µ + ∂ µ ϕ a /e to cancel the ϕ a term in Eq. (8) without changing the form of S ef f,a , resulting in L a also has an effective Lorentz invariance, which we also preserve against quantum corrections and take to be a symmetry of the effective action S ef f,a . We can further perform a chiral gauge transformation ψ a → e iΦaγ 5 /2 ψ a to cancel the phase of the Dirac mass in Eq. (14). Owing to the valley U (1) gauge invariance and the effective Lorentz invariance, Fujikawa's method suggests that the Jacobian of the measure (Eq. (13)) would contain a topologically nontrivial factor [71,72], which enters into the effective action as where F a,µν = ∂ µ A a,ν − ∂ ν A a,µ , and "..." includes all other terms. In this work, we only consider the leadingorder linear response to A and u zz . Through an explicit evaluation of Feynman diagrams in Appendix. C, we find that the only leading-order linear response contained in "..." is the trivial correction to the permittivity and permeability in the material, which can be absorbed into the Maxwell term of A. Hence, all nontrivial leading-order linear responses come from the first term of Eq. (15). After omitting all of the higher-order and trivial terms, we can split Eq. (15) into three parts by using Eq. (10): S ef f,a = S ef f,a,θ + S ef f,a,Σ + S ef f,a,uzz .
S ef f,a,θ is the action for valley-separated axion electrodynamics where the valley axion field is given by the phase of the CDW order parameter as S ef f,a,Σ takes the form of a Chern-Simons theory and describes a valley layered QAHE with a valley Hall conductivity given by the CDW wavevector Q Because the two valleys are related to each other by TR symmetry, then Σ H,a (a = 1, 2) take opposite values in each of the two valleys and add to a net-zero total Hall conductivity, coinciding with the vanishing Chern number required by the TR symmetry. In general, an odd-integer valley layered QAHE indicates that a gapped Weyl-CDW with well-defined valleys is in a WTI phase if all higher energy bands are topologically trivial, as discussed and numerically confirmed in Sec. IV. The last term S ef f,a,uzz describes the strain-induced effect and reads where B y = F 31 . The effective action (Eq. (16)) is one central result of this work. We emphasize that the validity of Eq. (16) replies on the fact that the system must be gapped everywhere -including the boundary -if the system is finite-sized [22].

C. Piezoelectric Effect and DPME
Because S ef f,a,θ and S ef f,a,Σ in Eq. (16) have opposite signs in each of the two valleys, then the contributions of S ef f,a,θ and S ef f,a,Σ sum to zero in the total effective action. Hence, the total effective action only includes S ef f,a,uzz , which can be rewritten as Eq. (22) indicates that the total action relies on a nonzero electron-strain coupling, implying that the response of the action characterizes the deviation of the electron from the homogeneous deformation of the sample [17].
The total current derived from S ef f can be decomposed into two parts j P E is the total low-energy valley Hall current induced by the pseudo-electric field which, as required by m y symmetry, lies in the xz plane. Eq. (11) and Eq. (20) together imply that the pseudoelectric field and the low-energy Hall conductivity both have opposite signs in the two valleys, such that the induced Hall currents add constructively to give a total nonzero value. (See Fig. 2(b).) In Eq. (24), j P E is thus the bulk-uniform piezoelectric current induced by the CDW wavevector Q, where the piezoelectric coefficient is given by Hence, j P E can be understood as a 3D stack of 2D valley Hall systems in which each layer exhibits the 2D piezoelectric effect discussed in previous literature [18][19][20][21]. We would like to emphasize that Eq. (25) only includes the low-energy contribution to the piezoelectric current, while the high-energy contribution to the piezoelectric current is typically also present in realistic materials. Nevertheless, the high-energy contributions to the piezoelectric current should also be uniform in the bulk of the system. In Eq. (23), j M takes the form of a magnetization current in which the total orbital magnetization M is induced by a pseudo-electric field through the valley axion field Physically, Eq. (26) can be understood from the bulkboundary correspondence as follows. First, given a gapped and symmetry-preserving boundary, the CDW phase φ smoothly changes from a constant value φ 0 in the bulk to zero in the vacuum, implying that the magnetization current j M is localized on the boundary. According to the bulk-boundary correspondence of the axion field, the surface valley Hall conductance (along the normal direction of the surface) should take the form θ bulk a 2π on any surface, where θ bulk a = (−1) a−1 φ 0 . Hence, the surface-localized magnetization current j M is simply the surface Hall current induced by the pseudoelectric field, as shown in Fig. 2(c). The surface current generates a uniform bulk magnetization of the form which is the DPME proposed in this work, as illustrated in Fig. 1(d). Unlike the piezoelectric current in Eq. (24) originating from the 2D valley Hall conductance, j m in Eq. (28) originates from the fundamentally 3D bulk valley axion field, and is fundamentally different from Eq. (24).

IV. MINIMAL TB MODEL FOR TR-INVARIANT WSM
The analysis in Sec. III is based on low-energy effective field theory. It is natural to ask whether our low-energy prediction of a DPME in TR-invariant Weyl-CDWs remains valid in the presence of high-energy bands (or equivalently in a UV completion). Furthermore, the lowenergy analysis in Sec. III relies on "valley" quantum numbers; however, valley index is neither a generic symmetry of TB models, nor a symmetry of real solid-state materials. To address these questions, we will construct a minimal TB model of a TR-invariant Weyl-CDW and compute the bulk-average value of the valley axion field and the DPME, which we will compare to those predicted by the effective action.
Prior to the onset of a CDW, we begin with an orthorhombic lattice, in which we choose for simplicity the lattice constants to be a x = a y = a z = a 0 . We then consider there to be two sublattices in each unit cell. We next place a Kramers pair of spinful s orbitals on one sublattice and a Kramers pair of spinful p y orbitals on the other sublattice, resulting in a four-component basis c † k,i,s , where i = 1, 2 is the sublattice index and s = ± is the spin index. We then construct a four-band TB model H T B,0 that preserves TR and m y symmetries, which has the following form where in which τ 0,x,y,z and σ 0,x,y,z have been redefined to respectively act on the sublattice and spin indices, and where explicit expressions for d 1,2,3,4,5 are provided in Appendix. D. Throughout this work, we will take the inverse of the lattice constant without strain 1/a 0 as the unit of energy, which occurs because we have employed a convention in which = c = 1. Using parameters specified in Appendix. D, H T B,0 hosts four Weyl points at that are related by TR and m y symmetries. The projection of the bulk Weyl points onto the (001) surface is shown with blue and red dots in Fig. 3(a), where a (red) blue dot indicates a Weyl point with a chiral charge (−1) 1. The topological surface Fermi arcs that connect Weyl points with opposite chiralities appear as bright curves in Fig. 3(a), and the linear dispersion of the (001)-projecting bulk bands from the Weyl points is shown in Fig. 3 On the k x = 0 plane, H T B,0 is gapped in the bulk and exhibits a nontrivial TR-protected 2D Z 2 index. The nontrivial Z 2 topology is indicated by the appearance of gapless helical modes along k x = 0 on the (001) surface, as shown in Fig. 3(c). We next add a CDW term that preserves TR and m y symmetries into the TB model, where the CDW coupling takes the form in which M 1 and M 2 are Hermitian matrices whose explicit forms are provided in Appendix. D, and where µ 1 and µ 2 are real scalar parameters. Eq. (32) suggests that the CDW term contains two channels that are characterized by two real coupling constants, µ 1 and µ 2 . Throughout this work, we will set |µ 1 + iµ 2 | = 0.3/a 0 for all of the numerical calculations for the TB model in the presence of the CDW. Unlike Ref. [73], we do not study the microscopic origin of the CDW order parameter in this work, as the main goal of introducing the TB model is simply to provide a UV completion of the low-energy theory on which our analysis is rigorously based. Owing to the lattice-commensurate nature of the CDW in H T B,CDW , H T B,0 + H T B,CDW has reduced (but not fully relaxed) lattice translation symmetry, where the new lattice constants of the modulated cell are given by a x = 2a 0 , a y = a 0 , and a z = a 0 . The CDW backfolds two Weyl points of the same valley index onto the same momentum in the reduced 1BZ to form an unstable 3D Dirac fermion, which then becomes gapped. The gap induced by the CDW is reflected by the appearance of a bulk gap at k x a x = π in Fig. 3(d), which stands in contrast to the gapless (WSM) bulk in Fig. 3(b).

A. Weak Z2 Topological Insulator Phase
While the bulk of the Weyl-CDW phase is gapped, Fig. 3(d) demonstrates the existence of the gapless helical edge modes along k x a x = π on the (001) surface. As shown in Fig. 3(e), there are also gapless helical edge modes on the (001) surface at k x = 0. Hence, the k x = 0 and k x a x = π planes both exhibit nontrivial TRprotected Z 2 topology. This indicates that the model H T B,0 + H T B,CDW is in a WTI phase characterized by a nontrivial weak Z 2 index vector (ν x , ν y , ν z ) = (1, 0, 0), where ν x , ν y , ν z are the weak Z 2 indices in the k x = 0, k y = 0, and k z = 0 planes, respectively.
The nontrivial WTI index vector can be understood from the odd-integer valley layered QAHE S ef f,a,Σ in the effective action Eq. (16), provided that the high-energy bands are trivial. To see this, we first project the TB model into each of the two valleys. From this, we see that H T B,0 + H T B,CDW reproduces Eq. (3) and Eq. (6) with and where ϕ is the U (1) gauge degree of freedom of the eigenvectors (further discussed in Appendix. D). We emphasize that, at this stage, we have not yet incorporated the effects of dynamical strain, which will be added in the next section. After restoring the Fermi velocities for Eq. (20) (Appendix. C), we can use Eq. (34) and Q = (π/a 0 , 0, 0) = (2π/a x , 0, 0) to derive the valley Hall conductivity induced by the CDW wavevector, which we find to be given by Eq. (36) implies that the two valleys as a set contribute two counterpropagating chiral edge modes for each k xindexed plane in the 1BZ. The two couterpropagating chiral modes exhibit a TR-protected crossing at k x = 0 and k x = π/a x , indicating the presence of a nontrivial weak Z 2 index ν x = 1. We pause to compare the WTI phase of H T B,0 + H T B,CDW (Eq. (29) and Eq. (32)) to the WTI Dirac-CDW phase in Ref. [74]. In Ref. [74], which was revised to include TR-symmetric semimetal-CDWs during the final stages of preparing this work, the authors study a TR-invariant Dirac semimetal that is gappd by a CDW. Specifically, in Ref. [74], two 3D Dirac points become coupled by the CDW order parameter, and the gapped Dirac-CDW additionally respects spatial inversion symmetry when the phase of the CDW order parameter φ 0 = 0, π. The authors of Ref. [74] find that φ 0 = 0, π correspond to two distinct WTI phases (with nontrivial weak Z 2 indices) that differ by a fractional lattice translation in the modulated cell. Although the model in Ref. [74] appears to be similar to H T B,0 + H T B,CDW with the two valleys moved to the same momentum, the m y symmetry enforced in this work is essentially different from the inversion symmetry at φ 0 = 0, π in Ref. [74]. This can be seen by recognizing that m y in this work is a symmetry of the Weyl-CDW for all values of φ 0 , including φ 0 = 0, π. We do note that H T B,0 + H T B,CDW has an effective spinless-m x symmetry at arg(µ 1 +iµ 2 ) = 0, π (see Appendix. D). However, the spinless m x symmetry is an artifact of our simple TB model, which only occurs because all of the Weyl points are located at the k z = 0 plane. Because the effective action Eq. (8) allows the four Weyl points to move away from the k z = 0 plane -thus breaking any form of m x symmetry (effective or physical), the artificial spinless m x symmetry of our TB model does not affect any of the conclusions of and analysis performed in this work.

B. Valley Axion Field
In addition to the valley QAH term S ef f,a,Σ , the term S ef f,a,θ in Eq. (16) indicates that the bulk average value of the valley axion field is after restoring the Fermi velocities (Appendix. C). To compare Eq. (37) with the TB model, we first set θ bulk a = 0 at arg(µ 1 + iµ 2 ) = 0 as the reference for the TB model.
We are required to set a reference for the evaluation of θ bulk a , because the valley axion field is defined over an open manifold, such that its bulk average value is reference-dependent. Hence, we may only evaluate the change of θ bulk a relative to a reference value. Because the low-energy result Eq. (37) suggests that θ bulk a = 0 for φ 0 = 0, we further set the U(1) degree of freedom ϕ in Eq. (35) to ϕ = 0 in order to match the zero points of θ bulk a in the low-energy result and in the TB model. We may then use φ 0 = arg(µ 1 + iµ 2 ) for the TB model and evaluate θ bulk a with [22] θ bulk where k is only integrated over the half of the reduced 1BZ that contains the ath valley. In Eq. (38), is the non-Abelian Berry curvature of the non-Abelian can be equivalently evaluated based on the Chern-Simons 3-form [22], this method would require to choose the gauge and k-space boundary condition very carefully, and thereby we use gauge-invariant Eq. (38) in this work.
The resulting numerical computation of θ bulk a for the a = 1 valley is shown in Fig. 3(f), and extremely closely matches the value θ bulk 1 = −φ 0 expected from the lowenergy action. However, there is still a quantitative deviation between the low-energy and TB results, which occurs because the bulk valley axion field is not defined over a closed manifold, and is thus not quantized, implying that high-energy degrees of freedom (which are necessarily present in a solid-state material) can drive the value away from the low-energy result. Nevertheless, the relatively small deviation between the TB and low-energy results in Fig. 3(f) suggests that the effect of high-energy modes on the valley axion field is small in the TB model employed in this work.
C. Piezoelectric Effect and the DPME Finally, we use the TB model to verify the straininduced piezoelectric effect and the DPME described by the term S ef f,a,uzz in Eq. (16). We incorporate the effects of strain into the TB model by adding a prefactor of (1 − (∆rz) 2 (∆r) 2 u zz ) for each hopping term in H T B,0 , where ∆r is the displacement of the hopping [75] (see Appendix. B and Appendix. D for further details). As a result, the low-energy projection of the extra strain term implies that the added strain reproduces Eq. (7) with Because the TB model is a WTI with the weak Z 2 index vector (1, 0, 0), we consider a slab configuration with N layers perpendicular to x with periodic boundary conditions along y and z. The slab can be viewed as a quasi-2D system, and we may therefore calculate the 2D piezoelectric tensor of the slab using [17] where A slab n,i = −i ϕ k y ,k z ,n |∂ k i |ϕ k y ,k z ,n and |ϕ k ,n is the periodic part of the Bloch states of the slab Hamiltonian in the presence of strain. By inserting a projection operator onto each layer of the slab, we can then derive the 2D piezoelectric tensor for each layer, which we label .., N is the layer index (see Appendix. D for details). The 2D current density induced by the infinitesimal dynamical strain for each layer is then given by Using this method, we next calculate the 2D z-directional current density of each layer for varying values of φ 0 and N = 20 by setting ϕ = 0 in Eq. (35) and using φ 0 = arg(µ 1 + iµ 2 ). We note that the current along the y direction conversely vanishes in our numerics at each value of φ 0 , due to the bulk m y symmetry at each value of φ 0 .
In Fig. 4(a), we plot the current density distribution As schematically shown in Fig. 4(b), we can decompose the current density distribution into a uniform background current j 2D z (averaged over the layer index) and a layer- The uniform background current j 2D z characterizes the uniform piezoelectric response, and, as shown in Fig. 4(c), the piezoelectric current is nearly independent of φ 0 , as expected from the low-energy expression Eq. (26).
On the other hand, the layer-dependent contribution to the layer current density δj 2D z (l x ) is asymmetrically distributed. Specifically, δj 2D z (l x ) exhibits opposite signs near the two surfaces, resulting in a bulk magnetization M bulk y . To calculate the bulk magnetization, we treat the l x th layer as a uniform 2D system that is infinite in the y, z directions but finite in the x direction as . From this, we then express the 3D current density as where Θ l x (x) = 1 for x ∈ [(l x −1)a x , l x a x ] and Θ l x (x) = 0 otherwise. We next take the magnetization at the center of the sample x 0 = a x N/2 derived from the Biot-Savart law to be the bulk magnetization, which yields where the uniform background contribution naturally vanishes. In Eq. (43), we have chosen sgn(0) = 0, be- where M 0 = eu zz /a 0 . As shown in Fig. 4(d), the TB and low-energy results are of the same order of magnitude (the deviation is smaller than 70% of the low-energy result). This agrees with our earlier determination that the high-energy modes have relatively small effects on the valley axion field (Fig. 3(f)). In particular, the TB and low-energy results in Fig. 4(c) match extremely well as φ 0 approaches ±π (the deviation is smaller than 7% of the low-energy result). As discussed below, the agreement between the TB and low-energy results can be attributed to the TB model exhibiting boundary gap closings at exactly φ 0 = ±π.

V. BOUNDARY TQPT AND DPME JUMP
In this section, we will first show that the slab configuration of H T B,0 + H T B,CDW has a boundary gap closing at φ 0 = ±π, which we will show to be a boundary TQPT that changes the surface Z 2 index and induces a discontinuous change of the DPME. This boundary gap closing accidentally happens on the two surfaces of the slab at the same value of arg(µ 1 + iµ 2 ), allowing us to fully interpret the discontinuous change of DPME within the low-energy theory. We will then add an extra term to H T B,0 + H T B,CDW to split the accidental simultaneous surface gap closing, resulting in a more realistic model in which the gap closings on the two surfaces occur at different values of arg(µ 1 + iµ 2 ). Lastly, we will demonstrate that a jump in the DPME still occurs across each surface gap closing, though the low-energy theory is incapable in fully describing the jump due to the unavoidable presence of a gapless boundary helical mode on one side of the jump. Throughout this section, we will continue to choose φ 0 = arg(µ 1 + iµ 2 ) by setting ϕ = 0 in (35), except in Sec. V C.

A. TB Model
According to Eq. (35), φ 0 only appears in the TB model as cos(φ 0 ) and sin(φ 0 ) in µ 1 and µ 2 , respectively, and thus any TB result must be periodic in φ 0 . Hence, tuning µ 2 from negative to positive while keeping µ 1 < 0 should drive φ 0 from −π to π and give a jump of the magnetization, as shown in Fig. 4(d). The dramatic difference between the current distributions at φ 0 = ±0.9 in Fig. 4(a) provides evidence of the expected jump in the bulk magnetization. Fig. 5(a) suggests that the jump of (a) the DPME at φ 0 = ±π happens along with the boundary gap closing while the bulk stays gapped. Moreover, the gap closing manifests as one 2D gapless Dirac cone in each valley on each surface perpendicular to x, as shown in Fig. 5(b-d). Because there are two TR-related Dirac cones on one surface at the gap closing, then the surface gap closing has the same form as the 2D Z 2 transition that happens at a TR-related pair of generic momenta [67,76,77]. Because the bulk remains gapped across the transition, then the surface gap closings represent examples of boundary TQPTs, which can be detected by jumps in the DPME. Another signature of the boundary TQPT appears in the domain wall structure shown in Fig. 6(a). We consider a slab configuration that is split into two parts along the z direction, where each part exhibits a different value of the CDW phase φ 0 (specifically φ + 0 for z > 0 and φ − 0 for z < 0), while all other parameters in the slab are taken to be the same for the two parts. The sample is set to be periodic in the y direction and open in the x direction. In our numerical calculations, we have specifically employed a slab with 20 layers along x and have fixed φ + 0 = 0.8π.
We note that we have also chosen 20 layers along z for z < 0 and 20 layers along z for z > 0, but we do not depict the additional surface modes at large |z| (i.e. the leftmost and rightmost surfaces in Fig. 6(a)).
We plot the phase diagram of the domain wall structure in Fig. 6(b) by varying φ − 0 from φ + 0 to π and then from −π back to φ + 0 . As a result, we identify two phases, which we label as I and II. The phase I contains the point φ − 0 = φ + 0 , implying that the surfaces of both parts of the slab are not related by boundary TQPTs, such that boundary between the two domains is gapped. As φ − 0 is varied from π + 0 − to −π + 0 + , the gap closes on the top and bottom surfaces normal to the x direction in the z < 0 region, as discussed above and shown in Fig. 5. This indicates that the z < 0 region has undergone a pair of boundary TQPTs to enter phase II. The appearance of 1D gapless helical modes at the edge of the z = 0 interface in phase II confirms the presence of a nonzero surface relative Z 2 index for the two sides of the gap closing (Fig. 6(c)). The 1D gapless helical modes in phase II persist until φ − 0 reaches φ + 0 − π, where the gap closes in the 2D bulk of the interface as shown in Fig. 6(d). The interface gap closing again manifests as two gapless Dirac cones at two valleys and thus changes the Z 2 index of the interface, coinciding the disappearance of the helical edge modes in phase I.

B. Low-energy Effective Theory
We now interpret the boundary TQPT in the TB model from the perspective of the low-energy theory. According to Eq. (8), one bulk Dirac cone has two mass terms, and thus the bulk gap closing for Eq. (8) requires fine-tuning at least two parameters, which typically does not occur in a realistic model or material. On the (100) surface, the projections of the valleys are along the m y Tinvariant line, and the gap closing along this line only requires fine-tuning one parameter, according to Ref. [21]. The analysis in Ref. [21] further suggests that the gap closing appears as one gapless surface Dirac cone for each valley (Fig. 7(a)) and is thus a surface Z 2 transition, coinciding with the TB results in Fig. 5 and Fig. 6(c). The parameter values for which the gap closings appear depend on the boundary condition that we choose in Eq. (8) (see Appendix. C for a special boundary condition that realizes both surface gap closings at φ 0 = π). Nevertheless, the codimension-1 nature of the gap closing indicates that, even if the boundary conditions are varied, it is still difficult to remove the gap closing point. When the boundary conditions are changed, the gap closing instead shifts to a different value of φ 0 . Indeed, as we will shortly show in using a TB model with an extra term that splits the simultaneous boundary gap closing, the boundary phase transitions are movable in arg(µ 1 + iµ 2 ), but globally unremovable. This agrees with the picture presented in Ref. [74], in which tuning φ 0 pumps 2D TI layers in the WTI phase until a layer reaches the system boundary, causing a surface gap closing. The same argument can also be applied to the (100) surface.
In general, the gap closings on the (100) and (100) surfaces do not happen at the same critical value of φ 0 . We find that simultaneous surface gap closings only occur when the system configurations (parameter values or boundary conditions) are designed in a fine-tuned manner such that unrealistic (i.e. artificial) effective symmetries appear in the effective action (such as an effective TR symmetry within one valley after omitting the CDW wavevector). This suggests that the presence of simultaneous surface gap closings in the above TB model H T B,0 + H T B,CDW is accidental. In this accidental (finetuned) case, a simultaneous gap closing changes the surface valley Hall conductance by ±e 2 /2π, which, according to the bulk-boundary correspondence of the axion term, results in a change of the θ bulk a -or equivalently φ 0 -by 2π. Combined with Eq. (28), the 2π jump of φ 0 further results in a jump of the magnetization  gapless Dirac cones at the interface of the domain when the two sides of the domain wall differ by π in φ0. In general, the number of 2D gapless Dirac cones at the interface is 2 + 4n where n is a non-negative integer. Fig. 4(d). Therefore, the boundary TQPT and the induced jump of the DPME can be captured within the low-energy theory when valleys are well-defined and when the gap closings happen simultaneously on both surfaces. Lastly, we will use the low-energy theory to explain the gap closing in the 2D bulk of the interface of the domain wall (Fig. 6(d)). Within the low-energy theory, if two gapped Dirac cones have a mass phase difference π and form a domain wall structure, then there must be an odd number of 2D gapless Dirac cones localized at the interface [73]. Therefore, when φ ± 0 in Fig. 6(a) differ by π, an odd number of 2D gapless Dirac cones appear at the interface for each valley, as shown in Fig. 7(b), and the gap closing correspondingly changes the Z 2 index of the interface. Indeed, the description given by the low-energy theory coincides with the TB result shown in Fig. 6(b) and (d).

C. Separate Gap Closings on Two Surfaces
In the final part of this section, we will discuss the more general (and also more realistic) case in which the gap closings on the (100) and (100) surfaces occur at different values of arg(µ 1 + iµ 2 ). In this subsection, we do not set ϕ = 0 in Eq. (35) and thus in general φ 0 = arg(µ 1 + iµ 2 ). We can shift the gap closings on two surfaces in opposite directions in arg(µ 1 + iµ 2 ) by adding an extra TR-and m y -symmetry-preserving term H extra in the TB model, as schematically shown in Fig. 8(a). (The explicit form of H extra is provided in Appendix. D). Because each surface gap closing is a Z 2 TQPT, then the Z 2 index of the entire slab is changed across the transition, resulting in two phases with different Z 2 indices for the whole slab (i.e. the slab as a whole, for varying arg(µ 1 +iµ 2 ), is or is not a 2D Z 2 TI). Exactly which phase of the slab has nontrivial Z 2 index is determined by the number of layers, owing to the nontrivial weak Z 2 index in the bulk (Fig. 8(a)). To numerically model this more generic case, we choose appropriate parameter values for the model (Appendix. D) to split the simultaneous boundary phase transition into two boundary transitions: one at arg(µ 1 + iµ 2 ) = 0.9π on the (100) surface and another at arg(µ 1 + iµ 2 ) = −0.9π on the (100) surface. This results in the appearance of two phases in the slab: one phase that includes arg(µ 1 + iµ 2 ) = 0, and another that includes arg(µ 1 + iµ 2 ) = π. As shown through slab Wilsonloop calculations [74] in Fig. 8(b,d), the 5-layer slab is Z 2 trivial in the arg(µ 1 + iµ 2 ) = π phase, and is nontrivial in the arg(µ 1 + iµ 2 ) = 0 phase. On the other hand, Fig. 8(c,e) indicate that the 6-layer slab is Z 2 trivial in the arg(µ 1 + iµ 2 ) = 0 phase and nontrivial in the arg(µ 1 + iµ 2 ) = π phase. The 5-layer and 6-layer results can be generalized for all odd-layer and even-layer slabs, respectively, as long as the number of layers is large enough to avoid any additional layer-dependent gap closings. In general, this is consistent with the early recognition of odd-even boundary modes in WTIs [80], and with the picture established in Ref. [74] in which a TRinvariant Weyl-CDW phase can be captured by a stack of 2D TIs whose normal vectors lie parallel to the wavevector Q, where the position of the 2D TI in each cell is set by arg(µ 1 + iµ 2 ).
The DPME predicted by the effective action is valid only when the slab is Z 2 trivial, as the nontrivial Z 2 index of the slab necessarily indicates the presence of a gapless helical mode on the side surface, violating the gapped boundary requirement (i.e., the validity of DPME predicted by the low-energy effective action requires a gapped and symmetry-preserving boundary). The failure of the effective action can also be seen from the bulk-boundary correspondence. The bulk-boundary correspondence of axion electrodynamics implies that an unambiguous bulk value of the axion field (in the units of 2π) should be equal to the Hall conductance of every gapped surface (in unit of e 2 /(2π)) [22]. When one surface undergoes a Z 2 transition, the valley Hall conductance on that surface changes by e 2 /(2π), while the valley Hall conductance on the other surface remains constant. As a result, at least on one side of each surface transition, different surfaces infer different bulk values of the valley axion field, indicating the incapability of the effective action in predicting the DPME [22]. The underlying physical reason for the failure of the effective action is that there is a contribution from the gapless helical mode to the DPME that cannot be captured by the lowenergy effective action. Therefore, in general, it is not always appropriate to set ϕ = 0 in Eq. (35). Instead, one should choose ϕ such that the DPME predicted from the low-energy action matches the TB result when the slab is Z 2 -trivial.
Nevertheless, the discontinuous change of the DMPE should still exist across the gap closing on one surface; the jump will just contain two contributions in the more realistic case. One contribution arises from the appearance of an extra gapless helical mode. The other contribution is given by the discontinuous change of the surface valley Hall conductance (or more directly, the discontinuous change of the strain-induced surface current). The slab configuration in Fig. 8(a) allows us to demonstrate the second contribution to the DPME jump ( Fig. 9(a)), because the periodic boundary conditions along y and z avoid the contribution from the side-surface helical modes. As long as the effective action is valid in the slab-Z 2 -trivial region, the total change of the DPME over the slab-Z 2 -nontrivial region can still be predicted. Moreover, in a domain wall structure like Fig. 6(a), the gapless surface domain wall mode still appears across the gap closing on one surface of the z < 0 side, where the domain wall mode is exactly the extra gapless helical mode that corresponds to the change in the slab Z 2 index on the z < 0 side of the domain wall.

VI. CONCLUSION AND DISCUSSION
Using low-energy theory and TB calculations, we have in this work introduced the DPME of TR-invariant WSMs in the presence of a bulk-constant (static and homogeneous in the bulk) CDW that gaps the bulk Weyl points. The DPME is a fundamentally 3D strain effect that specifically originates from a valley axion field. We further demonstrate a discontinuous change of the DPME across a boundary Z 2 TQPT by tuning the phase of the CDW order parameter. The discontinuous change of the DPME can serve as a bulk experimental signature of the boundary TQPT. Although we have only considered a pair of TR-related valleys, the analysis performed in this work can straightforwardly be generalized to multiple pairs of TR-related Weyl points, as long as one does not enforce additional crystalline symmetries that restrict the total DPME to be zero.
To probe the DPME, one can measure the induced magnetic field outside of a Weyl-CDW sample. However, the uniform piezoelectric current also generates a magnetic field outside of the sample, and thus it is important to devise a means of distinguishing the uniform piece from the DPME. One solution is to measure the magnetic field just outside of the (010) surface of a sample with a (100)-directed CDW that gaps the bulk Weyl points. As shown in Fig. 9(b), the magnetic field given by the DPME in this geometry is directed along y, whereas the magnetic field induced by the uniform piezoelectric current is directed along x. Additionally, the helical modes on the surfaces perpendicular to y and z can be gapped out by finite-size effects (similar as the gapped side surfaces of (a) The strain-induced bulk magnetization calculated from the realistic slab configuration of HT B,0 + HT B,CDW + Hextra (see Fig. 8). The number of layers is chosen to be 20. The two black dashed lines indicate the two surface transitions at arg(µ1 + iµ2) = ±0.9π. (b) A sample with finite size along three directions. Outside of the (010) surface, the uniform piezoelectric current generates a magnetic field BP E along x, while the magnetic field BDP M E given by the DPME is directed along y. The CDW wavevector is oriented along the x direction. the axion insulators [47] and antiferromagnetic topological insulators [37,38]), when the phase of the CDW order parameter and the number of layers are chosen to guarantee that the 2D Z 2 index of the slab is trivial. The helical modes can alternatively be removed by a side-surface CDW [81]. Based on Eq. (44), we estimate the order of magnitude of the response coefficient for the DPME to be |∂M /∂u zz | ∼ 0.8e/Å for ξ y ∼ 1eV [20], φ 0 ∼ π, and v y ∼ 10 −4 c [82]. We find that the DPME response coefficient has the same units as the 2D piezoelectric coefficient, which is reasonable, because the 3D bulk magnetization in the DPME manifests as the strain-induced 2D surface current density. Hence, we may directly compare the value of the DPME response coefficient with typical experimental values of 2D piezoelectric coefficients (∼ 10 −20 C/Å) [83], suggesting that the above estimated value of the DPME response coefficient is experimentally observable. An important direction of future study is to formulate the contribution from the helical mode to the DPME, which is relevant in the case in which there is one gapless helical mode left on the side-surface for certain values of the phase of the CDW order parameter with respect to the number of layers (see Fig. 8).
In addition to a direct probe of the DPME, our theory predicts the existence of a boundary TQPT and the appearance of 1D gapless helical modes along surface domain walls of the CDW phase. The gapless helical domain wall fermions can in principle be probed through scanning tunneling microscopy. Moreover, it is intriguing to ask whether the boundary TQPT separates two 3D phases with different boundary-obstructed topology [84], and to elucidate the precise relationship between the boundary TQPT and symmetry-enhanced topological surface anomalies [66,81,85].
In this work, we have focused on TR-invariant gapped Weyl-CDWs. TR-invariant WSMs have been realized in a number of non-centrosymmetric systems, including NbAs [86], TaAs [87][88][89][90], and (TaSe 4 ) 2 I [82,91], and an axionic Weyl-CDW phase has recently been demonstrated in (TaSe 4 ) 2 I [51]. We emphasize that the intuitive picture of the DPME given in Sec. II is applicable as long as the low-energy physics of a given system is well-captured by 3D Dirac fermions with complex mass terms. This implies that the DPME may also exist in other 3D Dirac materials [92]. In the current work, we have treated the axion field (or the CDW phase) as a fixed background field. After taking into account the dynamics of CDW (such as a phason), the effective action in Eq. (21) suggests a nonzero coupling between the phason of the CDW and the strain field, implying the intriguing possibility of strain engineering the CDW phase angle in TR-invariant WSMs. In this section, we derive the mean-field CDW term in Eq. (6). Here, we use the imaginary time and allow the temperature to be nonzero.
We first convert Eq. (3) to the imaginary time, resulting in where q = (ω, q), ω = (2n + 1)π/β is the fermionic Matsubara frequency, dω = (2π/β) ω if temperate T is not (2π) 4 e ixq+ir·ka,α c q,a,α , and xq = ωτ + q · r. We consider two channels of the interaction and where g 1 > 0 and g 2 < 0, and g 2 is a perturbation with |g 2 | ∼ 0. The g 1 term is just a double copy of that used in Ref. [49]. As shown below, g 1 accounts for the nonzero CDW magnitudes while g 2 determines their relative phase.
To derive the mean-field CDW term, we first perform the Hubbard-Stratonovich transformation on S int,1 Then, we have (A5) Now, we perform the mean-field approximation. We neglect the quantum fluctuation of m a (x) and only consider a classical m a (x) that minimizes S M F [ m], meaning that m a (x) ∝ ψ † x,a,− ψ x,a,+ . Comparing m a (x) with Eq. (6) of the main text, we can define m a (x) = m a e i(−1) a−1 Q·r where m a is the CDW parameter mentioned in main text. Since we only care about the CDW order parameter that is constant in the bulk of the system, we choose m a to be independent of x. As a result, we arrive at a simplified Z: where , and V is the total volume of the system.
Next, we derive S M F to the first order of g 2 . The g 2 -independent part of S M F , labeled as S M F,0 , reads which gives Define m a = |m a |e iφa , then we see that S M F,0 does not depend on φ a . For the first order of g 2 , we have As a result, we have S M F = S M F,0 + S M F,1 to the first order of g 2 . Next, we minimize S M F . First, for φ a , we have Since g 2 < 0 and I a > 0, φ 1 + φ 2 = 2nπ minimizes S M F . The 2π ambiguity will disappear after introducing a gapped boundary, and a symmetry-preserving boundary would give in the bulk. Second, for |m a |, we have where ω 2 + i q 2 i v 2 i ≤ Λ 2 is used. The equation for |m a | with v x = v y = v z = 1 matches that in Ref. [49], which indicates that we need to have a large enough g 1 to have the nonzero CDW magnitude. The solution to the above equation has the form where |m 0 | is independent of g 2 . In the main text, we directly neglect the g 2 in |m a | and choose |m a | = |m 0 |. Then, where λ ∈ R labels the spin-orbit coupling, where [U g ] G s ,Gs = δ G ,RG e −iG ·t (R S ) s s and R S is the matrix representation of R in the spin subspace. Then, we can define g u = {R|t u = (1 + u)t} and have For the TR symmetry that acts on the Bloch states, we have With these symmetry representations, we can derive the symmetry properties of h u (k) from those of h 0 (k), which are useful for the change (ii). Suppose [g, H 0 ] = 0. Then, we have g u (H R −1 uR )g −1 u = H u based on Eq. (B1) and Eq. (B3). As a result, we have For TR symmetry, suppose [T , H 0 ] = 0, and then we have T H u(t) T −1 = H u(t) , which gives

Tight-binding Model
The above formalism in general is hard to deal with analytically. More commonly, we deal with the TB model, which reads in the absence of deformation, or reads with deformation. Here c † R+τ i is a vector whose components stand for orbital, spin, etc, and are labeled by β i as c † R+τ i,βi . Note that β i can take different ranges of values for different sublattices. In this work, we approximate where δ = R + τ i − R − τ i is treated as a column vector. Since the above expression only involves the symmetric part of u, this only takes in to account the effect of strain. It is reasonable sicne a global rotation of the system cannot induce any response. Let us define with N redefined as the total number of lattice sites. From {c † r,β , c r ,β } = δ rr δ ββ , we can derive that c † k,i,βi and c † k,i,βi have the same anti-commutation relation With Eq. (B20), the Hamiltonian can be re-expressed as Therefore, similar as the general formalism, the strain effect to the TB model includes (i) c † k → c † k and (ii) h 0 (k) → h u (k) , where c † k and c † k have the same commutation relation, and h u=0 (k) = h 0 (k). The similarity also exists for the symmetry properties. First, Eq. (B11) for the lattice translations still holds here. Second, if [g, H 0 ] = 0 for a space group operation g, then c † R+τ i furnishes a representation of g, i.e., with M M i i g is defined to be zero for i and i being different kinds of atoms. Third, if [T , H 0 ] = 0 for the TR operation T , then with M i T the representation of T in the β space, and Eq. (B14) and Eq. (B16) hold here for a different definition of U T : Therefore, the strain effect is formally the same for the general Hamiltonian and the TB model.

Low-Energy Model
In this part, we project the general Hamitonian or the TB model to the low energy subspace. We consider a group of orthonormal vectors v α (k a ) that satisfy where a, α are re-defined to label the valley and energies. We choose v α (k a )'s so that they furnish a representation of the symmetry group of H 0 : if g is a symmetry of H 0 , and where q only takes a small symmetric neighborhood around k a . As a result, b † q and b † q have the same commutation relation and the effective Hamiltonian reads Here we use the form for the general Hamiltonian, since the derivation for the TB model is equivalent. Clearly, the effect of deformation again includes with (U ef f g ) a α ,aα = [W a g ] α α δ k a ,Rka e −iRka·t if g is a symmetry of H 0 , and The symmetry properties of b q can be derived by limiting u → 0 in the above expression.
In general, h ef f 0 (q) and h ef f u (q) can be derived from h 0 (q) and h u (q) using the perturbation theory, respectively. However, this is not always straightforward to be done analytically, so sometimes we derive their form from symmetries.
Note that H ef f 0 and H ef f u should have the same symmetry properties as H 0 and H u , respectively. Then,we have Now we restore the original definition of a, α, and consider the case discussed in the main text. With this scheme and the following symmetry representations we can obtain the leading-order electron-strain coupling as Converting to the field operator c t,q,a,α and using ψ t,r,a,α = e ika,α·r d 3 q (2π) 3 e iq·r c t,q,a,α , we can obtain Eq. (7) of the main text. Here the low-energy approximation allows us to extend the range of q to R 3 and treat a, α as internal indices. Since the strain cannot change the commutation relation of the field operator, the measure of the functional integral in the partition function does not change with the strain.
Appendix C: More Details on the Low-Energy Theory In the section, we provide more details on the effective action and the boundary condition of the low-energy theory.

Effective Action
In this part, we derive Eq. (16) in the main text from Eq. (14) and Eq. (12) of the main text. Let us first simplify Eq. (14) of the main text. As mentioned in the main text, |m| = |m(r)| equals to a constant |m 0 | in the bulk, and then we can define M (r) such that (s.t.) |m(r)| = |m 0 | + M (r). Since M (r) is constantly zero in the bulk regardless of the gapped boundary, it cannot contribute to any bulk response, and thus we can neglect it in the following derivation, simplifying Eq. (14) of the main text to From the above equation, a generic term in the effective action has a specific nonnegative powers of A a , A a,5 , and Φ a , labeled as n 1 , n 2 , and n 3 , respectively. Among all these terms, the charge response comes from those with n 1 ≥ 1, since the functional derivative of the action with respect to A is required to derive the current. We next adopt the leading order approximation, which consists of three parts. First, since we only consider the currents that are linear in u zz or A, the terms in the effective action that we are interested in must have n 1 ≤ 2. Second, we only keep terms to the linear order of φ and Q for simplicity, and then we can simplify Eq. (C1) to The first two approximations further require n 2 ≤ 2 and n 3 ≤ 1 and forbid (n 1 , n 2 , n 3 ) = (1, 2, 1), (2, 1, 1), (2, 2, 0), (2, 2, 1). As a result, we have only eight possible values of (n 1 , n 2 , n 3 ) that might have nonzero contribution to the effective action within the first two approximations, whose Feynman diagrams are summarized in Fig. 10. The third approximation is that since u zz and A are chosen to slowly vary along with t and/or r, we keep at most two space-time derivatives of them.
In the Feynman diagrams, the fermion propagator reads 1 i S(k) with since this form leaves the measure invariant. The Fourier transformation of A a , A a,5 , and Φ a follows the same rule as above. In this section, we adopt the real time: x µ = (t, r) µ and define q µ = (ν, q). In the following, we evaluate each graph. For the derivation, we do not use the dimensional regularization due to the existence of Levi-Civita symbol, but adopt the classic Adler's method for the chiral anomaly, which does not choose a regularization scheme as discussed in Ref. [71 and 72]. Eq. (C2) has effective Lorentz invariance and U(1) gauge invariance, and we preserve it to all orders of quantum correction.
Before evaluating each diagram, we would like to discuss a subtlety. The U (1) gauge field is added as ∂ x µ → ∂ x µ +ieA µ (x) for ψ t,r . According to Eq. (B41), ψ t,r actually corresponds to a fermion at x phy = (t, r phy ) = (t, (1+u)r). Therefore, A µ (x) is related to the actual physical U (1) gauge field A phy µ (x phy ) as Nevertheless, we can directly replace A(x) by the physical U(1) gauge field A phy (x) in this work as discussed below. All current responses can be split into two classes depending on whether the response involves the electron-strain coupling parameter ξ i or not. For the response that involves ξ i , it must at least involve the electron-strain coupling ξ i u zz of power one, which means the strain effect in A(x) (if appears) would make the response non-linear in u zz . Since we at most consider the response to the first order of u zz , A(x) should be directly replaced by A phy (x) for the response that involves ξ i . The ξ i -independent response would stay unchanged even if the electron-strain coupling limits to zero. If we keep the strain effect in A(x) for this type of response, a ξ i -independent strain-induced current might appear. Such a strain-induced response corresponds to the motions of the electrons that exactly follow the homogeneous deformation, which is ambiguous at the linear order for the infinitely large system according to Ref. [17] and thus must be neglected. Therefore, we should also directly replace A(x) by A phy (x) for the response that does not involve ξ i . In sum, we can treat A(x) as A phy (x) in our work. In other words, all the strain-induced linear current responses derived here characterize to what degree the electrons fail to follow the homogeneous deformation, which is what we should consider according to Ref. [17].
However, there is a tricky part here. Eq. (C6) is just one way to assign momentum to the graph Fig. 10(a), and there are infinite many other ways as with p independent of k. As Tr [γ µ S(k)] k 2 →∞ − −−− → 1/k and we do not choose any regularization scheme, we have d 4 k Tr [γ µ S(k)] has U V divergence d 4 k1/k. As a result, for nonzero p, meaning that the expression for Fig. 10(a) is ambiguous [71,72]. This ambiguity only appears for the UV divergence faster than logrimathiric divergence, which also appears for some other graphs in Fig. 10 as discussed below.
Nevertheless, we can use symmetry and physics to remove this ambiguity for Fig. 10(a). Eq. (C6), if nonzero, breaks the effective Lorentz invariance and suggests that a nonzero current would exist without any external perturbation. Therefore, it should be restricted to zero, meaning that Eq. (C6) is the only allowed way of assigning momentum.
Then, for f ρ (k) where we use which can be derived with Wick rotation. Here g µν = diag(−1, 1, 1, 1) µν is the metric. Owing to Eq. (C15), the expression of R ρ (q) in Eq. (C12) should be after taking into account the ambiguity induced by the momentum shift. Nevertheless, we can evaluate R ρ up to a momentum shift. To do so, we need to introduce the widely used trick called Feynman parametrization [71]: With this trick, R ρ becomes 2 , p ρ marks the ambiguity given by the shift of the momentum, including the change from k to l in the integral also contributes the momentum shift. The effective Lorentz in-variance requires R ρ (Λq) = Λ ρ ρ R ρ (q) for Lorentz transformation Λ. This symmetry requirement imposes a constraint on the momentum shift of the R ρ (q), p ρ must have the form p ρ = B 1 (q 2 )q ρ , which partially removes the ambiguity. As a result, we have Substituting the above expression into Eq. (C12), we find Eq. (C12) is zero.
Using Eq. (C11), the above expression can be evaluated to zero.
d. Aa Aa Term Fig. 10(d) formally is the same as the loop correction to the photon propagator for quantum electrodynamics [71], which has the form where C 0 is a constant. According to the expression of F a , summing the above expression over a gives where the terms that cannot contribute to the charge response have been neglected. Therefore, this term stands for the trivial correction of the permittivity and permeability inside the material. and At large k 2 , we have and Tr where we use Tr [γ µ1 γ µ2 ...γ µn ] = 0 for odd n .
f. AaAa,5Aa,5 Terms Fig. 10(g) reads where U µνρ q1,q2 = 1 2 and the shift of k in the integral of f µνρ U (k, q 1 , q 2 ) causes ambiguity as i 24π 2 (p ν g µρ + p ρ g µν + p µ g νρ ) . (C37) Here we use and Then, we can evaluate the integral of f µνρ U (k, p 1 , p 2 , p 3 ) up to a shift of integration variable: and the Taylor expansion with respect to the p i is used for the last equality. Substitute the p 1 = q 1 , p 2 = 0, and p 3 = −q 2 in the above equation and use we arrive at where the leading order term would vanish after symmetrizing by q 1 ↔ q 2 and ν ↔ ρ. We can neglect the O(q 3 ) terms in the above equation since the leading order approximation allows at most two space-time derivatives on the fields, meaning that we only need to consider the terms up to q 2 order. Then, since Eq. (C37) indicates that Eq. (C34) only holds up to a shift in integration variable, we have where p µ labels the total ambiguity, and p µ can depend on q 1 or q 2 but must be invariant under q 1 ↔ q 2 . This ambiguity is removed by the U (1) gauge invariance (Ward identity) that requires which further requires that p should be set to zero. As a result, we conclude that Fig. 10 where V µνρ q1,q2 = i 2 and Since f µνρ V (k, p 1 , p 2 , p 3 ) decays as 1/k 3 for large k 2 as f µνρ V (k, p 1 , p 2 , p 3 ) its integral is ambiguous as where Eq. (C11) and Eq. (C38) are used. It means Eq. (C47) holds up to a shift of momentum. Now we evaluate the integral of f µνρ V (k, p 1 , p 2 , p 3 ): where the first equality uses Eq. (C18) and Eq. (C41), the second equality uses Eq. (C11) and the Tyler expansions with respect to p i , the third equality uses Eq. (C42) and and the last equality uses Eq. (C50). Substituting p 1 = q 1 , p 2 = 0 and p 3 = −q 2 into the above equation and neglect O(q 3 ) terms according to the leading order approximation, we can derive the expression of V from Eq. (C47)as where p δ should depend on q 1 , q 2 and changes its sign under q 1 ↔ q 2 . Then, in general p δ = q 1,δ C 1 (q 1 , q 2 ) + q 2,δ C 2 (q 1 , q 2 ), where C 2 (q 2 , q 1 ) = −C 1 (q 1 , q 2 ). Again, the ambiguity is removed by the U (1) gauge invariance which requires C 1 = C 2 = 0. Eventually, we know Fig. 10(h) is negligible within the leading order approximation.
Since f µν T (k, p 1 , p 2 , p 3 ) the integral of f T has at most logarithmic divergence, and thus is unambiguous under the shift of integration variable. Then, we can evaluate the integral where the first equality uses Eq. (C18) and Eq. (C41), the second equality uses Eq. (C11) and the Tyler expansion with respect to p i , and the third equality uses dF n 1 = 1. Substituting p 1 = q 1 , p 2 = 0, and p 3 = −q 2 and neglecting the O(q 3 ) terms, we have which means the leading-order contribution from Fig. 10(i) to the effective action has the form of the axion term: Only this term contains nontrivial contribution to the leading-order linear response.

i. Restoring Fermi Velocities
As shown above, only Fig. 10(d) and (i) have nonzero contribution to the leading-order linear response. Fig. 10(d) only gives the correction to the permittivity and permeability in the material and thus is trivial. The only nontrivial leading-order linear response comes from Fig. 10(i), which gives Eq. (16).
Note that the following transformation is used to derive Eq. (8) and Eq. (14) of the main text, Then, we can perform the inverse transformation on Eq. (16) to restore the Fermi velocities. As a result, Eq. (11) of the main text becomes Eq. (18) of the main text becomes Eq. (20) of the main text becomes Eq. (21) of the main text becomes S ef f,a,uzz = sgn(v x v y v z ) eξ y 4π 2 v y dtd 3 r(φ + Q · r)u zz B y , Eq. (25) of the main text becomes and Eq. (28)  The derivation here suggests that the axion term comes from i|m 0 |Φ aψa γ 5 ψ a term, while the A a,5,µψa γ µ γ 5ψ has zero contribution. This seems to contradict the fact that A a,5,µψa γ µ γ 5 ψ a in the |m 0 | = 0 case should account for the axion term of the chiral anomaly. The reason for this seeming contradiction is that the Taylor expansion with respect to the momentum of A a /A a,5 /Φ a that we perform above is invalid in the |m 0 | → 0 limit. In another word, the method used here does not have the proper |m 0 | → 0 limit and thus cannot restore the massless chiral anomaly. However, we still adopt this method since it can reproduce the previous experimentally-verified results [49,51] and the results derived above are verified by the TB model as discussed in Appendix. D.
k. Momentum-Cutoff Correction to Eq. (C55) The nontrivial contribution to the leading order response is given by Eq. (C55). In the above, we have derived the response without imposing a finite momentum cutoff. In this part, we will impose a momentum cutoff |k| < Λ in Eq. (C55) to discuss the possible correction. Λ here is different from that in Appendix. A.
In the above derivation, we use Eq. (14), in which the momentum origin of the fermion field is typically not at the Weyl point. However, the order of magnitude of the derived correction to the response coefficient will be left invariant under a shift of the origin to the Weyl points, since such a shift can only change the forms of A a,5 and Φ a , and therefore cannot affect the form of the response coefficient according to Eq. (C73).

A Symmetry Preserving Boundary Condition
In this part, we present a symmetry-preserving boundary condition that can realize the boundary TQPT. We consider a slab sample with open boundary perpendicular to x and thickness 2x B , i.e., |x| < x B , while the momenta along y and z are kept as good quantum numbers. For the discussion of the boundary condition, we can omit the gauge field since it has no effect on the boundary condition. Moreover, we can also neglect strain since we choose it to be homogeneous in the entire space, meaning that the difference between the material and vacuum is solely accounted Owing to the monotonicity of Im(m), the φ0 should take values in (−π, π) and the continuous φ cannot take the path that passes ±π as |x| increases. (c) shows the spatial dependence of m for Im(m0) = 0 and Re(m0) < 0, where a gapless mode appears at the boundary.
for by the spatial dependence of the CDW order parameters. We choose the boundary to preserve the TR and mirror symmetries, and thereby we may focus on one valley since the other one is related by the TR symmetry. In this case, we can freely rotate the fermion bases to cancel the ϕ a and Q terms in Eq. and Im[m] to be smooth and monotonic for simplicity (see Fig. 11(a) for an example), and thus |m| is always continuous.
Owing to Im[m(|x| → ∞)] = 0, we can always set φ(|x| → ∞) = 0 for the vacuum. Furthermore, we always try to choose a continuous φ. When Im[m(|x| = 0)] = Im[m 0 ] = 0, the monotonicity requires Im[m] cannot take zero values for finite |x|. As a result, a continuous φ can only take values in (−π, π) and thus requires φ(|x| = 0) = arg(m 0 ) ≡ φ 0 to only take values in (−π, π), since φ otherwise must pass ±π as |x| increases and breaks the monotonicity of Im[m].(See Fig. 11(b).) When Im[m 0 ] = 0, the monotonicity requires Im[m(|x|)] = 0 for any x and φ must be a step function with φ(|x| < x B ) = nπ with n odd and φ(|x| > x B ) = 0. This discontinuity of φ for Im[m 0 ] = 0 comes from the gapless boundary mode given by the real mass domain wall Re[m 0 ] < 0 and Re[m(|x| → ∞)] → +∞, as well as the vanishing Im[m] = 0, as shown in Fig. 11(c). When φ 0 = π (e.g., Fig. 11(a)), the nonvanishing Im[m] guarantees the boundary to be gapped and makes sure that φ is continuous. Since the derivation of the response from the effective action is only valid when the system is gapped everywhere and φ is continuous, Eq. (28) of the main text is only valid for φ 0 = π, which gives us a uniform strain induced magnetization deep in the bulk of the system. If we keep Re[m 0 ] < 0 and tune Im[m 0 ] from 0 − to 0 + , φ 0 should jump from −π + 0 + to π + 0 − , leading to a jump of bulk magnetization. This jump is induced by the gap closing at the boundary of the system, while the bulk of the system stay gapped.
In the above discussion, we choose smooth monotonic functions for m. The existence of the gapless boundary mode for each valley and the magnetization jump is stable against any symmetry-preserving perturbation, as long as the two valleys are well defined. It is because the gapless mode for one valley on one surface is a Weyl point in the (q y , q z , φ 0 ) space, meaning that the perturbations can only shift the appearance of the gapless mode and the magnetization jump to other values of φ 0 instead of removing them. It coincides with the fact that the gap closing only needs 1 fine-tuning parameter.

Appendix D: Details on TB Model
In this part, we build a TB model to reproduce the results derived from the effective action.

Without CDW
We first consider the case without CDW. The model is built on a square lattice with lattice constant set to a 0 , and each lattice site contains two sub-lattice atoms, one at τ 1 = (0, 0, 0) and the other at τ 2 = (0, 0, 1/2)a 0 , as shown in Fig. 12(a). We put a spinful s orbital at τ 1 and a spinful p y orbital at τ 2 , and the bases read |R + τ i , s with s =↑↓ the spin index. According to the Fourier transformation the representations of the symmetry operations read where c † R+τ i,s is the creation operator for |R + τ i , s , and τ and σ label the sublattice and spin indices, respectively. With certain nearest-neighbor hopping terms, we choose the following symmetry-allowed form for the strained TB model: where the k dependence of d's is implied, the strain-induced redefinition of c k discussed in Appendix. B is implied, and d 1 = n 0 − 1 + cos(k x a 0 ) + n 2 cos(k y a 0 ) + (1 − u zz ) cos(k z a 0 ) where u zz stands for the normal strain along z. H T B,0 is just H T B,u with u zz = 0. Eigenenergies take the form ± d 2 1 + d 2 3 + ( d 2 2 + d 2 4 ± |d 5 |) 2 , and we consider half filling, resulting in the gapless condition d 1 = d 3 = d 2 2 + d 2 4 − |d 5 | = 0.