Orbital Edelstein effect from density-wave order

TL;DR

This paper demonstrates that a charge density wave order in a spin-orbit free electron model induces a temperature-dependent orbital Edelstein effect by reducing symmetry to a gyrotropic class.

Contribution

It introduces a model showing how density-wave order can generate a bulk orbital magnetization via the orbital Edelstein effect without spin-orbit coupling.

Findings

Orbital Edelstein effect appears in a spin-orbit free model with density-wave order.

The effect is temperature dependent and linked to symmetry breaking.

The model transitions into a gyrotropic crystal class enabling the effect.

Abstract

Coupling between charge and spin, and magnetoelectric effects more generally, have been an area of great interest for several years, with the sought-after ability to control magnetic degrees of freedom via charge currents serving as an impetus. The orbital Edelstein effect (OEE) is a kinetic magnetoelectric effect consisting of a bulk orbital magnetization induced by a charge current. It is the orbital analogue of the spin Edelstein effect in spin-orbit coupled materials, in which a charge current drives nonzero electron spin magnetization. The OEE has recently been investigated in the context of Weyl semimetals and Weyl metals. Motivated by these developments, we study a model of electrons without spin-orbit coupling which exhibits line nodes that get gapped out by via symmetry breaking due to an interaction-induced charge density wave order. This model is shown to exhibit a temperature dependent OEE, which appears due to symmetry reduction into a gyrotropic crystal class.

Equations

$${\mathbf{s}}_{\mathbf{k}n}=\frac{ge}{2m_{\text{e}}}⟨u_{\mathbf{k}n}|\mathbf{s}|u_{\mathbf{k}n}⟩,$$

$${\mathbf{m}}_{\mathbf{k}n}=\frac{e}{2\hslash }\text{Im}⟨{\nabla }_{\mathbf{k}}{u}_{\mathbf{k}n}|\times (H_{\mathbf{k}}-{\epsilon }_{\mathbf{k}n})|{\nabla }_{\mathbf{k}}{u}_{\mathbf{k}n}⟩,$$

$$\mathbf{M}=\frac{1}{\mathcal{V}}\sum _{\mathbf{k}n}{f}_{\mathbf{k}n}\left({\mathbf{m}}_{\mathbf{k}n}+{\mathbf{s}}_{\mathbf{k}n}\right),$$

$$f_{\mathbf{k}n}=f_{\mathbf{k}n}^0+e\tau (\mathbf{E}\cdot {\mathbf{v}}_{\mathbf{k}n}){\frac{\mathrm{d}f}{\mathrm{d}\xi }|}_{\xi ={\epsilon }_{\mathbf{k}n}-\mu },$$

$$M_{\kappa }={\alpha }_{\kappa \lambda }{E}_{\lambda },$$

$${\alpha }_{\kappa \lambda }$$

$$=\tau \frac{e}{\hslash }\sum _{\mathbf{k},n}{\frac{\mathrm{d}f}{\mathrm{d}\xi }|}_{\xi ={\epsilon }_{\mathbf{k}n}-\mu }(m_{\mathbf{k}n,\kappa }+s_{\mathbf{k}n,\kappa })v_{\mathbf{k}n,\lambda },$$

$$K=\sum _{\mathbf{k}}{{\psi }_{\mathbf{k}}}^{†}\left(d_{\mathbf{k}}^0+{\vec{d}}_{\mathbf{k}}^{}\cdot \vec{\tau }\right){\psi }_{\mathbf{k}},$$

$$d_{\mathbf{k}}^0$$

$$d_{\mathbf{k}}^1$$

$$d_{\mathbf{k}}^2$$

$$d_{\mathbf{k}}^3$$

$$V=\frac{1}{2}\sum _{(i,\alpha )\ne (j,\beta )}{V}_{ij}^{\alpha \beta }{n}_i^{\alpha }{n}_j^{\beta },$$

$${\mathbf{Q}}_0$$

$${\mathbf{Q}}_2$$

$${\rho }_{\mathbf{Q}}^{\text{s}}$$

$${\rho }_{\mathbf{Q}}^{\text{a}}$$

$$V_1$$

$$\alpha =\left(\begin{array}{ccc}{\alpha }_{aa}& & \\ & -{\alpha }_{aa}& \\ & & 0\end{array}\right).$$

$$\nu :={\tilde{\rho }}_{{\mathbf{Q}}_0}^{\text{a}},$$

$$d_{\mathbf{k}}^0$$

$$d_{\mathbf{k}}^1$$

$$d_{\mathbf{k}}^2$$

$$d_{\mathbf{k}}^3$$

$$d_{\mathbf{k}}^0$$

$$d_{\mathbf{k}}^1$$

$$d_{\mathbf{k}}^2$$

$$d_{\mathbf{k}}^3$$

$$d_{\mathbf{k}}^0$$

$$d_{\mathbf{k}}^1$$

$$d_{\mathbf{k}}^2$$

$$d_{\mathbf{k}}^3$$

$$m_{\mathbf{k}n,\kappa}=-\frac{e}{2\hbar}\varepsilon_{\kappa\lambda\zeta}\frac{\vec{d}_{\mathbf{k}}}{\big{|}\vec{d}_{\mathbf{k}}\big{|}^{2}}\cdot\left(\frac{\partial\vec{d}_{\mathbf{k}}}{\partial k_{\lambda}}\times\frac{\partial\vec{d}_{\mathbf{k}}}{\partial k_{\zeta}}\right),$$

$$d_{\mathbf{k}}^0$$

$$m_{\mathbf{k}n,a}$$

$$m_{\mathbf{k}n,b}$$

$$m_{\mathbf{k}n,c}$$

$$V_{\text{MF}}$$

$$=N\sum _{\mathbf{q}}^{\text{BZ}}\sum _{\alpha \beta }{V}_{\mathbf{q}}^{\alpha \beta }{{\rho }_{\mathbf{q}}^{\alpha }}^{\ast }{n}_{\mathbf{q}}^{\beta }.$$

Full text

I Introduction

The field of magnetoelectric effects has seen a revival in interest in the past decades on several fronts. The discovery and study of multiferroicity in correlated materials has uncovered unconventional mechanisms which can give rise to a large effective magnetoelectric coupling Spaldin and Fiebig (2005); Eerenstein et al. (2006); Ramesh and Spaldin (2007); Cheong and Mostovoy (2007); Spaldin et al. (2010); Lawes and Srinivasan (2011); Fuentes-Cobas et al. (2015); Chu et al. (2018). Similarly, the discovery of three-dimensional (3D) topological insulators has led to an exploration of novel magnetoelectric effects due to emergent axion electrodynamics in such topological phases Hasan and Kane (2010); Qi and Zhang (2011); Grushin and de Juan (2012); Pesin and MacDonald (2013); Schmeltzer and Saxena (2013); Mal’Shukov et al. (2013); Baasanjav et al. (2014); Morimoto et al. (2015); Xiao et al. (2018); Tokura et al. (2019). A part of the reason for the wide interest in such magnetoelectric effects partly stems from the technological potential of controlling charge degrees of freedom via magnetic fields or, conversely, tuning magnetic degrees of freedom via an applied electric field.

A prominent example of such a magnetoelectric effect is the nonequilibrium phenomenon of current-induced magnetization, which is also termed as kinetic magnetoelectric effect (KME)Levitov et al. (1985); Şahin et al. (2018). The intrinsic-spin variant of this effect, wherein a charge current in a spin-orbit-coupled conductor gives rise to bulk spin polarization and, hence, a net magnetization, is referred to as the Edelstein effect or the inverse spin-galvanic effectSinova et al. (2015); Manchon et al. (2015) and has been under study for several decadesGanichev et al. (2016). A great deal of experimental work has focused on the Edelstein effect in 2D systems, notably in thin-film semiconductorsKato et al. (2004); Silov et al. (2004); Kato et al. (2005); Sih et al. (2005); Stern et al. (2006); Yang et al. (2006) and at metal surfacesZhang et al. (2014). However, experiments on 3D materials have been scant, although some recent studies have reported its observation in trigonal tellurium Shalygin et al. (2012); Furukawa et al. (2017).

In recent years, it has come to light that 3D systems can have an intrinsic orbital contribution to the KME, analogous to the spin part and arising as a consequence of the orbital magnetic moment of Bloch bandsYoda et al. (2015); Zhong et al. (2016); Rou et al. (2017); Yoda et al. (2018); Şahin et al. (2018); Tsirkin et al. (2018); Flicker et al. (2018); Niu et al. (2018); Shi and Song (2019), notably in trigonal selenium and tellurium. Whereas the ordinary Edelstein effect (hereafter referred to as the spin Edelstein effect, SEE) relies on crystalline spin-orbit coupling (SOC) to give Bloch states a spin texture and, hence, is limited by the size of the SOC, the orbital Edelstein effect (OEE), also referred to as the inverse gyrotropic magnetic effectZhong et al. (2016), is determined solely by the geometry of the crystalYoda et al. (2015, 2018).

Chiral crystals are a subset of those that can exhibit the KME. Previous studies have considered trigonal selenium and tellurium and viewed their chiral nature as descending from a charge-density-wave (CDW) instability of a hypothetical parent phaseSilva et al. (2018), and others have studied optical gyrotropy as a probe for symmetry breaking in the chiral CDW phase of $1T$-$\mathrm{TiSe_{2}}$ Gradhand and van Wezel (2015) and in stripe-ordered cupratesOrenstein and Moore (2013). It has also been shown that Weyl nodes at the Fermi level can yield a large intrinsic contribution to the KME Yoda et al. (2015, 2018).

Our work builds on this theme, and explores the KME induced by symmetry breaking in a system with line nodes in the electronic band structure. The resulting phase is a non-chiral but gyrotropic crystal, and we study the concomitant temperature-dependent KME as a probe of the density-wave order. Below, we briefly review the OEE, before introducing our model Hamiltonian and presenting its theoretical study.

II Orbital Edelstein effect

Electrons in crystalline solids form Bloch bands with an intrinsic spin magnetic moment

[TABLE]

where $\mathbf{s}=\hbar\bm{\sigma}/2$ is the spin operator for the electrons. The modern theory of magnetization in solidsResta (2010); Thonhauser (2011); Vanderbilt (2018); Resta (2018) has discovered that such Bloch bands also host an intrinsic orbital magnetic moment given by

[TABLE]

where $H_{\mathbf{k}}$ is the Bloch Hamiltonian, $n$ is the band index, and $H_{\mathbf{k}}|u_{\mathbf{k}n}\rangle=\varepsilon_{\mathbf{k}n}|u_{\mathbf{k}n}\rangle$. This orbital magnetization has been shown to arise from the self-rotation of wavepackets in the semiclassical theory of electron dynamicsXiao et al. (2010); Thonhauser (2011). Since $\mathbf{m}_{\mathbf{k}n}\to-\mathbf{m}_{-\mathbf{k}n}$ under time reversal and $\mathbf{m}_{\mathbf{k}n}\to\mathbf{m}_{-\mathbf{k}n}$ under spatial inversion, it is clear that at least one of these symmetries must be broken in order for $\mathbf{m}_{\mathbf{k}n}$ to not be identically zero.

From the viewpoint of semiclassical dynamics, given an electron distribution function $f_{\mathbf{k}n}$, the instrinsic contribution to the net electronic magnetization is given byXiao et al. (2010); Resta (2010); Şahin et al. (2018)

[TABLE]

where $\mathcal{V}$ is the crystal volume. In thermodynamic equilibrium for a time-reversal symmetric system, the Fermi-Dirac distribution, $f_{\mathbf{k}n}^{0}=f(\varepsilon_{\mathbf{k}n}-\mu)$ forces zero a net magnetization $\mathbf{M}=0$ because of cancellation between contributions from opposite crystal momentaYoda et al. (2018). However, an asymmetric distribution function, such as that arising from an applied electric field, can generally give rise to nonzero net bulk magnetization.

Explicitly, to lowest order in an applied uniform DC electric field, the distribution function is modified asAshcroft and Mermin (1976)

[TABLE]

where $\mathbf{v}_{\mathbf{k}n}$ is the electronic group velocity, $\tau$ is the impurity-scattering relaxation time in relaxation-time approximation, and $e>0$ is the elementary charge. Hence, the magnetization arises as a linear response to an applied electric field,

[TABLE]

with the linear response tensor

[TABLE]

where $\kappa$ and $\lambda$ are Cartesian indices.

The form of the tensor $\alpha$ is significantly constrained by crystal symmetryYoda et al. (2018); Ganichev et al. (2016); Nye (1957); Landau and Lifshitz (1960). $\alpha$ is an axial rank-two tensor since it relates a polar vector, $\mathbf{E}$, to a axial vector, $\mathbf{M}$. Crystal classes whose point-group symmetries allow for nonzero axial rank-two response tensors are known as gyrotropic. The reason for this name is that the tensor governing natural optical activity, or gyrotropy, transforms in the same way as $\alpha$; thus KME and optical gyrotropy go hand in hand.

We note that the same symmetry constraints govern the appearance of nonzero spin and orbital contributions to $\alpha$, so both are expected to arise together, and there is no clear route to disentangling them in a 3D systemYoda et al. (2018); Furukawa et al. (2017). Indeed, the authors of Ref. Furukawa et al., 2017 conclude by speculating that the current-induced magnetization they observe in trigonal tellurium may be due not only to the well-known SEE, but also to the OEE.

III Model

As an example of the OEE brought about by symmetry breaking, we consider a tight-binding toy model of spinless fermions moving in a tetragonal crystal as shown in Fig. 1, consisting of identical atoms are arranged in layered square lattices. We assume a single isotropic orbital, and ignore the electron spin below; there are many cases where this is a useful starting point. In crystals with density-wave order driven by nearest-neighbor repulsion, such as we will consider below, both spin components behave in the same manner. Including spin then only leads to an extra factor-of-two in certain equations below. The other case where spin may be ignored is in spin-polarized systems which might be a useful description of states in a large energy interval around the Fermi energy in strong ferromagnets (i.e., in half-metals).

We define the nearest-neighbor (NN) lattice constant $a_{0}$, and lattice constant along the stacking-$\hat{\mathbf{z}}$ axis to be $c$. In the $x$-$y$ planes, we include NN hopping $t_{1}$ and next-nearest-neighbor (NNN) hopping $t_{1}^{\prime}$. In the $y$-$z$ and $z$-$x$ planes, we include NN hopping $t_{3}$ and the peculiar NNN hoppings depicted in Fig. 1, with $t_{2A}\neq t_{2B}$. This last choice differentiates staggered $A$ and $B$ sublattices, and endows the crystal with primitive translations $\mathbf{a}=a_{0}(\hat{\mathbf{x}}+\hat{\mathbf{y}})$, $\mathbf{b}=a_{0}(\hat{\mathbf{x}}-\hat{\mathbf{y}})$, and $\mathbf{c}=c\hat{\mathbf{z}}$. The crystal space group is $\mathit{P4/nbm}$, and its associated crystal class is $\mathit{4/mmm}$. Importantly, this crystal has centers of inversion at the middle point of every NN bond.

The inequivalence of the hopping amplitudes $t_{2A}$ and $t_{2B}$ for $\mathit{P4/nbm}$ symmetry may be rationalized as depicted in Fig. 1: if a secondary set of atoms is present only along the red bonds—which would be consistent with the space-group symmetry of the model—and their energy levels are far from the Fermi energy, their dynamics could be integrated out, with the end result of renormalizing the hopping amplitude between the primary atoms (shown in gray), effectively leading to $t_{2A}\neq t_{2B}$. The same reasoning holds if there are additional secondary atoms on the green bonds provided they are of a species different from those on the red bonds.

We use the following values for the hopping parameters: $t_{1}^{\prime}=0.7\,t_{1}$, $t_{2A}=0.1\,t_{1}$, $t_{2B}=0.4\,t_{1}$, and $t_{3}=0.5\,t_{1}$.

As is true for any crystal class with inversion symmetry, $\mathit{4/mmm}$ is non-gyrotropicNye (1957); Landau and Lifshitz (1960). However, if the symmetry of the crystal were to be reduced, for instance by the onset of CDW order, the inversion symmetry could be broken and, indeed, the ordered structure could fall into a gyrotropic class. Below, we first study the band structure of this model in the absence of interactions, before turning to the impact of nearest-neighbor repulsion.

III.1 Non-interacting band structure

The non-interacting Hamiltonian is $K=\sum_{i,j}-t_{ij}^{\alpha\beta}{c_{i}^{\alpha}}^{\dagger}c_{j}^{\,\beta}$, where ${c_{i}^{\alpha}}^{\dagger}$ creates an electron on sublattice $\alpha$ of unit cell $i$ (that is, on the atom at position $\mathbf{r}_{i}^{\alpha}$). In momentum space

[TABLE]

$c_{\mathbf{k}}^{\alpha}=N^{-1/2}\sum_{i}\mathrm{e}^{-\mathrm{i}\mathbf{k}\cdot\mathbf{r}_{i}^{\alpha}}c_{i}^{\alpha}$, $N$ is the number of unit cells, $\vec{\tau}$ is the vector of Pauli matrices (acting in sublattice space), and {\psi_{\mathbf{k}}}^{\dagger}=\Big{(}{c_{\mathbf{k}}^{A}}^{\dagger}\ \ {c_{\mathbf{k}}^{B}}^{\dagger}\Big{)}. Because of time-reversal symmetry, the hopping amplitudes $t_{ij}^{\alpha\beta}$ are necessarily real valued. We will measure momenta in units of inverse lattice spacing, and henceforth set $a=c=1$. With this, we arrive at

[TABLE]

The inversion symmetry $\mathit{\Pi}$ mentioned above combined with the Hamiltonian’s time-reversal symmetry $\mathit{\Theta}$ constrain $d^{3}_{\mathbf{k}}\equiv 0$ identically; for this reason, band touchings for this crystal will generically arise as line nodes, as seen in the non-interacting band structure (top panel of Fig. 6) Figure 2 shows the location of the line nodes in the Brillouin zone.

Spontaneous symmetry breaking (SSB), however, could change the crystal class to a less symmetric one. A simple scenario is a CDW phase in which the densities on the $A$ and $B$ atoms are inequal, in which case the space group becomes $\mathit{P\bar{4}2m}$, whose crystal class—$\mathit{\bar{4}2m}$—is gyrotropic.

III.2 Repulsive interactions

Next, we include repulsive interactions between NNs and next-nearest neighbors (NNNs); that is,

[TABLE]

where $n_{i}^{\alpha}={c_{i}^{\alpha}}^{\dagger}c_{i}^{\alpha}$ is the number operator for the atom $(i,\alpha)$ and we take $V_{ij}^{\alpha\beta}=V_{ji}^{\beta\alpha}$. We take nonzero repulsions $V_{ij}^{\alpha\beta}$ on the same bonds as the hopping amplitudes and use a similar naming scheme, as shown in Fig. 1. Note that, unlike for the hopping amplitudes $t_{2A}\neq t_{2B}$, we take the repulsion strengths to be $V_{2}$ on both the red and green bonds — we make this choice since it serves as a representative slice through the full parameter space and serves to illustrate some of the key ideas. Hence, the full Hamiltonian is given by $H=K+V$.

In the absence of spin order, a Hubbard term $U\sum_{i}n_{i\uparrow}n_{i\downarrow}$ would modify the MF Hamiltonian by a mere uniform offset in the chemical potential $\mu$. Hence, as we will restrict ourselves to ansätze without spin order, we do not include the Hubbard interaction.

We performed a MF calculation of the CDW order in the above system, allowing for a finite set of commensurate wavevectors. We identified the ordering wavevectors favoured by the interaction by considering a simple model of classical charges resting at each atomic site of the crystal of Fig. 1—see Appendix A. Based on this, we included in our ansatz the four ordering wavevectors

[TABLE]

The order parameters for the mean-field theory are the Fourier amplitudes $\rho_{\mathbf{Q}}^{\alpha}$ for the $\mathbf{Q}$ listed above, where $\alpha\in\{A,B\}$. However, we find it convenient to express the amplitudes on the $A$ and $B$ sublattices in terms of a symmetric and an antisymmetric part, respectively defined as

[TABLE]

Since $\rho_{\mathbf{Q}_{0}}^{\text{s}}=\rho$, the electron filling, we are left with 7 independent MFs.

IV Results

IV.1 Mean-field theory of charge-density-wave order

We focus on a cut of parameter space parametrized by $g$ such that

[TABLE]

According to the model of electrostatic charges, at these relative repulsion values (marked by a white star in Fig. 10), the interaction most favours an ordering with $\rho_{\mathbf{Q}_{0}}^{\text{a}}\neq 0$; however, the region of parameter space with $\rho_{\mathbf{Q}_{3}}^{\alpha}\neq 0$ most favoured is nearby (see Appendix A).

The MF calculation, for which the zero-temperature phase diagram is shown in Fig. 3, does indeed identify a swath of pure $\rho_{\mathbf{Q}_{0}}^{\text{a}}\neq 0$ order for $g$ sufficiently large and for approximately $0.3<\rho<0.6$—call this order phase I. In addition to this phase, we discover mixed phases at high and low filling in which both $\rho_{\mathbf{Q}_{0}}^{\text{a}}\neq 0$ and $\rho_{\mathbf{Q}_{3}}^{\alpha}\neq 0$, where $\alpha$ is either $A$ or $B$—call the order in these regions phase II.

For $g/t_{1}\gg 1$, we observe that the phase diagram is approximately symmetric under $\rho\rightarrow 0.5-\rho$; this is expected given the particle-hole symmetry of $V$, which becomes an approximate symmetry of $H$ when $K$ is nonzero but small compared with $V$.

The phase diagram at $\rho=0.5$ as a function of $T$ and $g$ is shown in Fig. 4 and reveals that raising the system temperature progressively supresses phase-I order.

IV.2 Magnetoelectric response in phase I

We studied the temperature-dependent SSB-induced OEE in broken-symmetry phase I, whose sole order parameter is $\rho_{\mathbf{Q}_{0}}^{\text{a}}$. Figure 5 shows the evolution of $\rho_{\mathbf{Q}_{0}}^{\text{a}}$ and the components of the response tensor $\alpha^{\text{orb}}$ at half filling ($\rho=0.5$) and fixed interaction strength ($g=1.5\,t_{1}$). The phase transition, which is seen to be continuous, takes place at $k_{\text{B}}T_{\text{c}}\approx 1.4\,t_{1}$.

As alluded to previously, the point-group symmetry in phase I is the $\mathit{\bar{4}2m}$, and a symmetry analysis reveals that symmetry constrains the response tensor $\alpha$ to the formNye (1957)

[TABLE]

As seen in the lower panel of Fig. 5, the calculated $\alpha^{\text{orb}}$ is indeed of this form. The nonmonotonicity of $\alpha^{\text{orb}}$ is understood via the MF band structure, seen in Fig. 6. The top panel shows the high-temperature ($k_{\text{B}}T=2.00\,t_{1}$) band structure. Although the system is in a metallic phase with many states within $k_{\text{B}}T$ (gray shading) from the chemical potential $\mu$ (dotted line), in the high-symmetry phase, $\mathbf{m}_{\mathbf{k}n}=0$ identically, so $\alpha^{\text{orb}}=0$. The middle panel shows the band structure at an intermediate temperature $k_{\text{B}}T=1.20\,t_{1}$ (near the peak in the $\alpha^{\text{orb}}$); the phase-I charge order has broken inversion symmetry and split the bands throughout the Brillouin zone, and the reduced symmetry of phase allows gyrotropic response, so $\alpha^{\text{orb}}\neq 0$. As the temperature is further lowered, the size of the band splitting increases and a full gap develops, with $\mu$ within this gap (lower panel of Fig. 5). Hence, it is clear that in the low-temperature limit, as the sum in Eq. 6 approaches a Fermi-surface integral, $\alpha^{\text{orb}}$ should again vanish.

IV.3 Role of line nodes

In this subsection, we inspect the contribution of the gapped-out line nodes to the magnetoelectric response tensor $\alpha^{\text{orb}}$ in phase I. As depicted in Fig. 2, in the high-symmetry phase, three independent line nodes exist in the first Brillouin zone: 1. the $X$-$R$ segment, 2. the $X$-$M$ segment, and 3. the $R$-$A$ segment—all other line nodes are related by symmetry.

In phase I, the charge order breaks the inversion symmetry $\mathit{\Pi}$; accordingly, in the MF Hamiltonian for phase I, $d_{\mathbf{k}}^{3}$ is no longer constrained to be zero, so the band touchings disappear. Rather, $d_{\mathbf{k}}^{3}$ is set by the phase-I order parameter $\rho_{\mathbf{Q}_{0}}^{\text{a}}$. As a shorthand, define

[TABLE]

where $\tilde{\rho}_{\mathbf{Q}_{0}}^{\text{a}}$, defined in Appendix B, is proportional to the phase-I order parameter $\rho_{\mathbf{Q}_{0}}^{\text{a}}$.

We expand about the (gapped) line nodes, denoting small deviations in momenum by $\tilde{\mathbf{k}}$. In the neighborhood of the (gapped) nodes, we have

$k_{a}=\mathrm{\pi}+\tilde{k}_{a}$, $k_{b}=\tilde{k}_{b}$, $-\mathrm{\pi}<k_{c}\leq\mathrm{\pi}$

[TABLE] 2. 2.

$-\mathrm{\pi}<k_{a}<\mathrm{\pi},\quad k_{b}=\mathrm{\pi}+\tilde{k}_{b},\quad k_{c}=\tilde{k}_{c}$

[TABLE] 3. 3.

$-\mathrm{\pi}<k_{a}<\mathrm{\pi},\quad k_{b}=\mathrm{\pi}+\tilde{k}_{b},\quad k_{c}=\mathrm{\pi}+\tilde{k}_{c}$

[TABLE]

We can compute the orbital magnetic moment about such a (gapped) line node: for a two-band model with Bloch Hamiltonian $H_{\mathbf{k}}=d_{\mathbf{k}}^{0}\tau^{0}+\vec{d}_{\mathbf{k}}\cdot\vec{\tau}$, $\mathbf{m}_{\mathbf{k}n}$ can be written in the closed formZhong et al. (2016); Yoda et al. (2018)

[TABLE]

where $\kappa$, $\lambda$, and $\zeta$ are Cartesian indices.

In the neighborhood of a (gapped) line node in the $\hat{\mathbf{a}}$ direction,

[TABLE]

where $f$, $g_{1}$, $g_{2}$, and $h$ are functions of $k_{a}$, and, hence, the components of $\mathbf{m}_{\mathbf{k}n}$ are given by

[TABLE]

where primes denote differentiation of the single-variable functions. To describe a line node in another direction, the indices in Eqs. 18 and 19 must be interchanged appropriately.

Figure 7 shows that the main contributions to $\alpha^{\text{orb}}$ originate from the the line nodes. The figure, for which we take $k_{\text{B}}T=1.20\,t_{1}$ as in the middle panel of Fig. 6, shows the orbital magnetic moment $m_{\mathbf{k}n,a}$ (left panel) as well as the product $-m_{\mathbf{k}n,a}v_{\mathbf{k}n,a}$ for the bottom ($n=1$, middle panel) and top ($n=2$, top panel) bands, which enter into the calculation of $\alpha^{\text{orb}}_{aa}$. [The component $\alpha^{\text{orb}}_{bb}$ is related by symmetry (Eq. 12), while all other components vanish.] These quantities are shown along slices $k_{c}=-\mathrm{\pi}$, $-\mathrm{\pi}/2$, [math], $\mathrm{\pi}/2$, and $\mathrm{\pi}$.

The form of $m_{\mathbf{k}n,a}$ in momentum space reflects all the symmetries of the point group $\mathit{\bar{4}2m}$ of phase I. The large values of $m_{\mathbf{k}n,a}$ are attributable to the presence of line nodes and are as expected fom Eqs. 19: for the slices $k_{c}=\pm\mathrm{\pi}/2$, we see the signature of the line node transverse to $\hat{\mathbf{a}}$, for which $m_{\mathbf{k}n,a}\rightarrow 0$ in proximity to the node, and for the slices $k_{c}=0$ and $k_{c}=\pm\mathrm{\pi}$, we also recognize contributions from the line nodes longitudinal to $\hat{\mathbf{a}}$, which are largest at the locus of the line nodes. The product $-m_{\mathbf{k}n,a}v_{\mathbf{k}n,a}$ still retains the general form of $m_{\mathbf{k}n,a}$, though with a modified symmetry that allows a nonzero integral over the Fermi surface.

The difference in the magnitude of the longitudinal contributions at $k_{c}=0$ and $k_{c}=\pm\mathrm{\pi}$ can be retraced to the different $g_{1}(k_{a})$ for the line nodes in those planes: while $g_{1}(k_{a})=2(t_{1}+t_{2A}+t_{2B})\cos(k_{a}/2)$ for the former, $g_{1}(k_{a})=2(t_{1}-t_{2A}-t_{2B})\cos(k_{a}/2)$ for the latter, and $|t_{1}+t_{2A}+t_{2B}|>|t_{1}-t_{2A}-t_{2B}|$ with our parameter values.

As shown by Eq. 6, up to constant factors, the integrand for $\alpha^{\text{orb}}_{aa}$ is $-m_{\mathbf{k}n,a}v_{\mathbf{k}n,a}$ multiplied by $-\mathrm{d}f/\mathrm{d}\xi|_{\xi=\varepsilon_{\mathbf{k}n}-\mu}$, a weighting factor concentrated within $k_{\text{B}}T$ around the chemical potential $\mu$. Hence, the equal-energy surfaces at $\varepsilon_{\mathbf{k}n}=\mu$ (solid black lines) and $\varepsilon_{\mathbf{k}n}=\mu\pm k_{\text{B}}T$ (dashed black lines) reveal the main contributions to the integral. At this temperature, for our choice of parameters, the largest contributions to $\alpha^{\text{orb}}$ are from the (gapped) line nodes at $k_{c}=\pm\mathrm{\pi}$.

V Discussion and conclusion

It is sometimes overlooked that gyrotropy can arise without chiral or time-reversal symmetry breakingOrenstein and Moore (2013); Ganichev et al. (2016). Since the point group $\mathit{\bar{4}2m}$ of the phase I order contains mirror symmetries (in our coordinates, one perpendicular to $\hat{\mathbf{x}}$ and another perpendicular to $\hat{\mathbf{y}}$), the structure of this broken-symmetry phase is an example of a non-chiral gyrotropic structure. This is unlike many of the examples of the OEE studied so farYoda et al. (2015, 2018), such as trigonal selenium and telluriumRou et al. (2017); Şahin et al. (2018). Furthermore, it constitutes a concrete example of longitudinal magnetization induced by a current in a mirror-symmetric structure, while it has been implied that this is not allowed by symmetry Yoda et al. (2015). While longitudinal magnetization is forbidden by mirrors perpendicular or parallel to the current, it is not forbidden by mirrors at $45^{\circ}$ angles from the current.

Similarly to a previously studied modelYoda et al. (2018), a picture of current flowing through solenoids provides a qualitative understanding of the OEE in phase I. Furthermore, this picture makes physically clearer why the OEE response vanishes (i) when $t_{2A}=t_{2B}$ — hence clarifying physically the necessity of assuming $t_{2A}\neq t_{2B}$ from the start —and (ii) in the absence of charge order, that is, when the atoms are indistinguishable. We imagine a current driven in the $\hat{\mathbf{b}}$ direction: the simplest solenoid-like paths that result in net displacement purely in the $\hat{\mathbf{b}}$ direction are shown in Fig. 8 and 9.

(i) Each row of atoms, whether on sublattice $A$ or $B$, is surrounded by a solenoid traced out of $t_{2A}$ (green) hoppings and another traced out of $t_{2B}$ (red) hoppings (Fig. 8). These two solenoids are of opposite helicites, and if $t_{2A}=t_{2B}$, their induced magnetic fields cancel identically, precluding any current-induced magnetizationwhether in the symmetric phase or in a charge-ordered phase. Indeed, if $t_{2A}=t_{2B}$, there exist additional mirror planes (perpendicular to $\hat{\mathbf{a}}$ and to $\hat{\mathbf{b}}$) that interchange the two solenoids.

(ii) If $t_{2A}\neq t_{2B}$, let us presume without loss of generality that $t_{2B}$ (in red) is dominant and ignore the $t_{2A}$ solenoids. Figure 9 shows that a solenoid whose central axis is a row of $A$ sites is of opposite helicity as one whose central axis is a row of $B$ sites. If the atoms on the two sublattices are indistinguishable, the magnetic fields induced by these two solenoids are equal and opposite and cancel when averaged over several lattice spacings. If the $A$ and $B$ atoms are distinguishable, however, the two sets of solenoids are distinct and can give rise to net magnetization. Indeed, while in the symmetric phase, there exist glide planes perpendicular to $\hat{\mathbf{a}}$ and along $\hat{\mathbf{b}}$ that interchange the dotted-line and solid-line solenoids.

It is surprising a priori that such a solenoid picture can exist in a nonchiral crystal (previous examples have focused on chiral structuresYoda et al. (2015, 2018); Şahin et al. (2018)); however, the two mirror symmetries of point group $\mathit{\bar{4}2m}$—one perpendicular to $\hat{\mathbf{x}}$ and the other perpendicular to $\hat{\mathbf{y}}$—do not bring the solenoids into themselves, but rather exchange the $a$-axis solenoids with the $b$-axis solenoids. The mirrors instead explain why the OEE coefficients along the $a$ and $b$ axis are opposite in sign, since the current is a polar vector and the magnetization is an axial vector.

In summary, we have discussed a simple model of a crystal with lines nodes, and shown how charge order can lead to a nonzero OEE. Such an OEE could be probed by nuclear magnetic resonance experiments. We have also discussed how the magnetoelectric response has a large contribution arising from the vicinity of line nodes. Our work suggests that lightly doped line-node semimetals in materials with weak spin-orbit coupling might be a promising place to search for large OEE.

Note

While completing this manuscript, we came across a recent theory preprint Ishitobi and Hattori (2019) which discusses the OEE induced by quadrupolar symmetry breaking in certain diamond lattice materials.

Acknowledgements

We thank D. A. Pesin for fruitful discussions. This research was funded by the Natural Sciences and Engineering Research Council of Canada and the Canadian Institute for Advanced Research. G. M. is supported by the Fonds de recherche du Québec - Nature et technologies. This research was enabled in part by support provided by Compute Ontario, Westgrid and Compute Canada. Computations were performed on the Niagara supercomputer at the SciNet HPC Consortium. SciNet is funded by: the Canada Foundation for Innovation; the Government of Ontario; Ontario Research Fund - Research Excellence; and the University of Toronto.

Appendix A Classical phase diagram

The ordering wavevectors favoured by the interaction of Eq. 9 were identified by determining the energy of CDW modes according to $V$ in a classical picture of electrostatic charges. In Fig. 10, we show which mode has the lowest energy (according to $V$) as a function of the (relative) sizes of the repulsion strengths $V_{1}$, $V_{1}^{\prime}$, $V_{2}$, and $V_{3}$. Note that, unlike in the main text, the modes are given in $x$-$y$-$z$ coordinates as \big{(}Q_{x},Q_{y},Q_{z}\big{)} (where we have set $a_{0}=c=1$), since it is not necessary to work with a multi-atom basis when considering only the repulsion term $V$. Furthermore, note that the yellow phase $(\mathrm{\pi},0,\mathrm{\pi})$ implicitly stands for itself as well as its symmetry-related counterpart $(0,\mathrm{\pi},\mathrm{\pi})$; the two are degenerate as expected by symmetry.

The ansatz for the mean-field theory described in the main text was chosen to potentially allow all the ground states that occur in this simplified model.

Appendix B Details of the mean-field calculation

Here, we provide additional information regarding the self-consistent MF calculation of the CDW order. The interaction of Eq. 9 was decomposed in the density channel, giving rise to the following MF interaction term:

[TABLE]

Here, $N$ is the number of unit cells in the crystal, $\rho_{\mathbf{q}}^{\alpha}=N^{-1}\sum_{i}\mathrm{e}^{-\mathrm{i}\mathbf{q}\cdot\mathbf{R}_{i}}\rho_{i}^{\alpha}$ (and likewise for $n_{\mathbf{q}}^{\alpha}$), and $V_{\mathbf{q}}^{\alpha\beta}$ is the Fourier transorm of $V_{ij}^{\alpha\beta}$, which is invariant under simultaneous translation of $i$ and $j$.

Choosing a closed set of commensurate wavevectors $\{\mathbf{Q}\}$ defines a reduced Brillouin zone (RBZ), which is mapped to the full Brillouin zone under addition of the wavevectors $\mathbf{Q}$. This property allows us to rewrite the integral of reciprocal-space-periodic functions as

[TABLE]

where the domain of the momentum sums is indicated above the summation symbol. Hence, we diagonalized the Hamiltonian $H_{\text{MF}}=K+V_{\text{MF}}$ in the RBZ by writing it in the form

[TABLE]

where $\tilde{\rho}_{\mathbf{Q}}^{\,\alpha}=\sum_{\beta}V_{\mathbf{Q}}^{\alpha\beta}\rho_{\mathbf{Q}}^{\,\beta}$ and $h_{\mathbf{k}}^{\alpha\beta}$ is the Bloch Hamiltonian, in our case given in Eq. 7: $h_{\mathbf{k}}=d_{\mathbf{k}}^{0}+\vec{d}_{\mathbf{k}}^{\vphantom{0}}\cdot\vec{\tau}$. Starting from a series of randomized values for the MFs $\rho_{\mathbf{Q}}^{\alpha}$, we iterated until the computed expectation values agree with the input MFs to within $10^{-6}$. We compared the Helmholtz free energy $F$ of the different ground states thusly obtained and selected the one with minimal $F$ at every point in parameter space.