Magnetic field induced symmetry breaking in nonequilibrium quantum networks

Symmetries in closed quantum systems have been a cornerstone in modern physics, introducing constraints to simplify and solve complex many-body systems. In recent years, inspired by their predictive capabilities, the basic principles of symmetries have been extended to open quantum systems. It has been proved that the dynamics of an open system can present several fixed points (asymptotic states) if the system is invariant under a symmetry [1] operation. This degeneracy of the dynamics generator determines the thermodynamic properties of the system and can potentially give rise to dynamical phase transitions [2, 3]. The symmetry related phenomena are surprisingly robust under weak symmetry breaking and are reflected in the metastable states [4, 5] instead of the asymptotic states.

Despite these novelties [6], little is known about how the breaking of symmetries affects the nonequilibrium properties of the system. In order to address this issue, we focus on cubic networks placed under a temperature bias, driving the system far from equilibrium. Cubic networks are one of the simplest closed systems with symmetries, found in many natural materials [7]. Recently, several experimental setups have been made accessible to study quantum transport in square lattices. Examples are optical lattices, where several transport experiments have been already performed [8, 9], and two-dimension ion traps with tunable spin–spin interactions [10–12]. Theoretically, due to the presence of open-system symmetries in square lattices both ballistic and diffusive subspaces can coexist [13, 14].

Inspired by these results and the simplicity of cubic networks, we aim to address if the closed-system symmetry can help us simplify the open-system nonequilibrium problem. We analyze the cubic networks driven far from equilibrium and obtain the underlying open-system symmetry operators and consequently the steady-state reduced density matrices that reflect the symmetries. Furthermore, we break the closed-system symmetry using Abelian magnetic fluxes [15] and show how this affects the open system leading to a control of nonequilibrium transport properties. In particular, we find that the direction of the magnetic field, and not its magnitude, plays a crucial role to control transport properties.

The rest of the paper is arranged as follows: in section 2 we describe the basic model of an open quantum cubic network in presence of a magnetic field. Section 3 reviews the basic notions of open-system symmetries and the non-triviality of obtaining the symmetry operators for open systems. We then treat analytically the simplest possible case of a system arranged as a square or cube in section 4. Here, due to the analytically obtainable symmetry operators, we get the initial condition independent steady states. We then in section 5 extend our model to quasi-one dimensional cubic networks and explore the transport numerically for different system-bath connections. In both the analytics and numerics we find that the magnetic field direction breaks symmetries that can dramatically affect the transport properties of the system. Finally, in section 6 we conclude.

The system of interest is a cubic network of non-interacting fermions subjected to an external magnetic field illustrated in figures 1–3, 6 and 7. The tight-binding representation of the system Hamiltonian is given by,

$$\begin{equation}H=-\sum _{r}{\varepsilon }_{r}{c}_{r}^{{\dagger}}{c}_{r}-\sum _{r,d=1,2,3}{t}_{d}{U}^{d}\left(\hat{r}\right){c}_{r+d}^{{\dagger}}{c}_{r}+\mathrm{h}.\mathrm{c}.,\end{equation} \tag{ 1 }$$

where r = (x, y, z) is the position of the lattice sites, ɛ_r is the onsite energy at position r, and t_d is the tunneling strength in the direction d, with d = 1, 2 and 3 corresponding to X, Y, and Z direction respectively. c^† and c are the fermionic creation and annhilation operators, respectively. The hopping element is isotropic. The onsite energies are different for the neighboring sites in the X–Y plane and homogeneous in the Z direction. This leads to the lattice having a period 2a (a being the norm of the primitive vector) in the X–Y plane causing the isotropic magnetic field to effect the couplings in all directions. The magnetic field effects are absorbed in the tunneling phases U¹(r) = 1, ${U}^{2}\left(r\right)=\mathrm{exp}\left[\text{i}2\pi B\left(x-y-\frac{1}{2}\right)\right]$, and ${U}^{3}\left(r\right)=\mathrm{exp}\left[-\text{i}\pi B\left(x-y\right)\right]$ with the cube length a = 1⁶. The phases are obtained under the Hasegawa gauge [16–20] that can be experimentally realized in artificial super-lattices [21–24] or cold atoms [15, 25, 26].

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Our system is connected to two reservoirs as shown in the illustrations of figures 1–3, 6 and 7 allowing us to explore the effects of magnetic fields on spatially symmetric nonequilibrium dissipative systems. We engineer the system-bath coupling Hamiltonian

$$\begin{equation}{H}_{S\alpha }=\sum _{l}{c}_{l}^{{\dagger}}{{\Pi}}_{i\ne l}{c}_{i}{c}_{i}^{{\dagger}}\sum _{k}{\nu }_{k}{d}_{k,\alpha }+\mathrm{h}.\mathrm{c}.\end{equation} \tag{ 2 }$$

with l index characterizing the sites connected to the reservoir α = L, R (L = left reservoir, R = right reservoir). It has been observed that engineering of baths and their coupling to ultracold atoms is possible and leads to interesting new features not observed in natural systems [27, 28]. The fermionic operators d, d^† correspond to the annihilation and creation operators of the reservoir, ${H}_{\alpha }={\sum }_{k}\,{{\epsilon}}_{k}{d}_{k,\alpha }^{{\dagger}}{d}_{k,\alpha }$, with ν_k being the system reservoir coupling strength. The specific form of the coupling allows the reservoirs to locally inject or extract a particle between the single particle states and the ground state. It is important to note here that for the extensively studied scenario of electron hopping between the system and reservoir ${H}_{S\alpha }={\sum }_{l}{c}_{l}^{{\dagger}}{\sum }_{k}\,{\nu }_{k}{d}_{k,\alpha }+\mathrm{h}.\mathrm{c}.$. This coupling Hamiltonian does not have any open-system symmetries and gives only one unique steady state. This is perhaps one of the crucial reasons that even for a simple cubic structure the effect of the magnetic field direction has not been observed in the literature.

The dynamics of the reduced density matrix (RDM) is governed by a Markovian dissipative Lindblad quantum master equation [29–32] that reads,

$$\begin{equation}\frac{\mathrm{d}\rho }{\mathrm{d}t}=\mathcal{L}\rho =-\mathrm{i}\left[H,\rho \right]+\sum _{\genfrac{}{}{0pt}{}{k=1,2}{\alpha =\text{L,}\;\text{R}}}{{\Gamma}}_{\alpha k}\mathcal{D}\left[{A}_{\alpha k}\right]\rho .\end{equation} \tag{ 3 }$$

Where $\mathcal{L}$ is the Liouvillian superoperator of the system, $\mathcal{D}\left[{A}_{\alpha k}\right]\rho ={A}_{\alpha k}\rho {A}_{\alpha k}^{{\dagger}}-\frac{1}{2}\left\{{A}_{\alpha k}^{{\dagger}}{A}_{\alpha k},\rho \right\}$, the Lindblad operators A_αk (${A}_{\alpha 1}={\sum }_{n}\,{c}_{n}^{{\dagger}}{{\Pi}}_{i\ne n}{c}_{i}{c}_{i}^{{\dagger}}\equiv {A}_{\alpha 2}^{{\dagger}}$) allow local injection (k = 1) or extraction (k = 2) of particles into the system. The sum over n in the Lindblad operators is restricted to n = 1, 5, 9, ⋯ (n = 1, ..., 4) for α = L and n = 3, 7, 11, ⋯ (n = N − 3, ..., N) for α = R for edge (face) connected baths. The injection rate Γ_α1 = Γf_α(₀) with ${f}_{\alpha }\left({{\epsilon}}_{0}\right)={\left(\mathrm{exp}\left[{\beta }_{\alpha }{{\epsilon}}_{0}\right]+1\right)}^{-1}$ being the Fermi–Dirac distribution. The ratio of the rates Γ_α1/Γ_α2 = exp[−β_α0] obeys local-detailed balance where β_α = 1/k_BT_α represents the inverse temperature of the αth bath and ₀ is the characteristic energy of the bath. The system-reservoir coupling ensures that without the magnetic field the generator can be split into invariant subspaces that lead to multiple steady-states [1, 3, 4, 6].

In this section we briefly review the relation between symmetry and invariant subspaces of the system dynamics (a more detailed explanation can be found in references [1, 3]). For a specific set of generators of the Markovian dynamics of our system $\left\{H,{\left\{{A}_{\alpha ,k}\right\}}_{\genfrac{}{}{0pt}{}{k=1,2}{\alpha =\text{L,}\;\text{R}}}\right\}$, if there is a unitary operator U s.t.

$$\begin{equation}\left[U,H\right]=\left[U,{A}_{\alpha ,k}\right]=0\quad \forall \enspace \left(\alpha ,k\right),\end{equation} \tag{ 4 }$$

then the system has a strong symmetry. In this case, we can prove that there are several invariant subspaces of the system dynamics, and therefore several orthogonal steady-states. To do so, we decompose the Hilbert space of the pure system, $\mathcal{H}$, by the spectral decomposition of U. If our symmetry operator has n_U eigenvalues we have that $U\vert {\psi }_{i}^{\left(k\right)}\rangle =\mathrm{exp}\left(\text{i}{\theta }_{i}\right)\vert {\psi }_{i}^{\left(k\right)}\rangle$ (i = 1, ..., n_U) and k = 1, ..., d_i with d_i being the dimension of the subspace corresponding to the eigenvalue exp(iθ_i). Using this we decompose the Hilbert space in the form $\mathcal{H}={\bigoplus}_{i}{\mathcal{H}}_{i}$, with ${\mathcal{H}}_{i}=$ span$\left\{\vert {\psi }_{i}^{\left(k\right)}\rangle :k=1,\dots ,{d}_{i}\right\}$. The operator space can be decomposed in a similar way as $\mathcal{B}\left(\mathcal{H}\right)={\bigoplus}_{ij}\,{B}_{ij}$, being B_ij = span $\left\{\vert {\psi }_{i}^{\left(k\right)}\rangle \langle {\psi }_{j}^{\left(l\right)}\vert :k=1,\dots ,{d}_{i};\right.$$\left.\,l=1,\dots ,{d}_{j}\right\}$.

We define now the left and right symmetry superoperators acting on the operators space as U_LO = UO and U_RO = O U, for any $O\in \mathcal{B}\left(\mathcal{H}\right)$. It is clear that the subspaces B_ij are the eigenspaces of these superoperators, meaning that U_LB_ij ⊆ B_ij as well as U_RB_ij ⊆ B_ij. As the system dynamics is generated by a superoperator $\mathcal{L}$ in the form of equation (3), and all the dynamics generators commute with the symmetry operator U it follows that the subspaces B_ij are invariant under the system dynamics $\mathcal{L}{B}_{ij}\subseteq {B}_{ij}$.

There is also the possibility of a superoperator S that commutes with the system Liouvilian $\left[S,\mathcal{L}\right]=0$ without having a strong symmetry associated. This is called a weak symmetry and we can show, by the same reasoning as with the strong symmetries, that a weak symmetry also implies the existence of invariant subspaces.

As we are working only with bounded systems we can apply Evans' theorem [39] to prove that in each subspace in the form B_ii there is at least one fixed point of the dynamics (eigenvectors of $\mathcal{L}$ with zero eigenvalue). It has also been proposed in reference [1] the possibility of finding fixed points in the subspaces ${\mathcal{B}}_{ij}$ with i ≠ j. These subspaces do not contain physical states as they include only zero trace operators. In this work we provide the first example of these non-physical fixed points.

The problem of finding the symmetry decomposition of a given system is still open. Thus, finding all the symmetry operators and invariant subspaces of an arbitrary system is a challenging task. In other words, without knowledge of all symmetry operators it is impossible to find all the invariant subspaces and consequently finding the RDMs within these subspaces is highly non trivial. Moreover, any naive attempt at numerically computing the steady-state RDMs will result in a matrix that is a linear combination of the true RDMs, and is generally not a true density matrix (see appendix B for more details on how to overcome this problem).

Given the cubic network in a magnetic field, we first provide the steady state solutions (ρ_ss = lim_t→∞  ρ) utilizing the underlying symmetries of the system. In the cases discussed here, using the symmetry operators we were able to obtain all the far from equilibrium steady states (states that are independent of initial condition), which is nontrivial to obtain in degenerated open quantum systems (see appendix B).

Figure 1 shows the various steady states for a single plaquette (square) illustrated in the first column. The edge coupled (top) and face coupled (bottom) reservoirs possess different open-system symmetries, giving rise to different steady states. The closed-system possess 4 symmetries (rotation by 0, π and flip about both diagonals) given by the plane symmetries of a square with different adjacent vertices. Not all of these four closed-system symmetries translate to open-system symmetries and the coupling of the system to the bath plays an important role (see equation (4)). In case of edge coupled plaquette, the open-system symmetry that survives in absence of magnetic field is the flip symmetry along the diagonal (the axis of symmetry is indicated by the dashed line in the illustration). The corresponding unitary symmetry operator is $U=\mathrm{exp}\left(\vert 2\rangle \langle 4\vert +\vert 4\rangle \langle 2\vert \right)$. The symmetry operator block diagonalizes the Liouvillian giving rise to two steady states as indicated in panels (a)–(d), obtained analytically and represented in the single particle basis. Out of the two states, one is a nonequilibrium steady state (NESS) that depends on the properties of the reservoir (panels (a) and (b)) and the other one is a pure state independent of the reservoir properties (panels (c) and (d)). The pure state turns out to be anti-symmetric under the exchange of states |2⟩ and |4⟩ as seen in panel (c). Such states carry zero current [33] and are also referred to as dark or compact localized states [34]. In presence of an anisotropic (z-direction) or isotropic magnetic field (panels (e) and (f)) the edge coupled system displays a single NESS that differs from the one obtained in the absence of magnetic field (panels (a) and (b)).

The face coupled single plaquette in the absence of magnetic field (panels (g)–(l)) gives three physical steady states with two pure states that are anti-symmetric with respect to exchange of sites |1⟩ − |3⟩ (panel (i)) and |2⟩ − |4⟩ (panel (k)). The corresponding symmetry operators are ${U}_{1}=\mathrm{exp}\left(\vert 2\rangle \langle 4\vert +\vert 4\rangle \langle 2\vert \right)$ and ${U}_{2}=\mathrm{exp}\left(\vert 1\rangle \langle 3\vert +\vert 3\rangle \langle 1\vert \right)$. In the absence of a magnetic field all steady states are real. Once the magnetic field (anisotropic or isotropic) is applied we break the flip symmetry along the diagonal but the π rotation symmetry still persists. Additionally, a new non-topological symmetry emerges that gives us three steady state solutions even in the presence of a magnetic field (panels (m)–(r)). In this case all the states are physical and we still find only one NESS. The two pure states (panels (o)–(r)) are non trivial, i.e., either symmetric or anti-symmetric w.r.t. the sites. Interestingly in this case the NESS is invariant under the influence of a magnetic field (in stark contrast to the edge coupled case).

The simplest single plaquette provides a range of unexpected results, like the presence of extra symmetries with the magnetic field, invariance of the steady state w.r.t. magnetic field, etc., which help describe the complete nonequilibrium picture. We have also obtained the RDMs for the cubic system (both edge and face coupled, see figures 2 and 3). The next example beyond a single plaquette, i.e., cube, increases the complexity drastically. The edge coupled ladder cube (cube without links between the plaquettes in the z direction) possesses 6 steady states for $\to {B}=\left(0,0,0\right)$, 5 steady states for $\to {B}=\left(0,0,B\right)$ (anisotropic field), and 1 for $\to {B}=\left(B,B,B\right)$ (isotropic field). The presence of an isotropic field destroys all open-system symmetries and gives a unique NESS, i.e., steady state that depends on the bath parameters. All the other states are pure states and are eigenstates of the system Hamiltonian.

Next we study the face connected cube in figure 3. In this case the direction of the magnetic field plays a significant role with 5 steady states being present without a magnetic field (panels (a)–(j)), 5 steady states for an anisotropic field $\to {B}=\left(0,0,B\right)$ (panels (k)–(t)) and only 1 state for an isotropic field $\to {B}=\left(B,B,B\right)$ (panels (u) and (v)). In each case, we have only one NESS and the rest of the states are pure states. Here even though the number of steady states without and with anisotropic field are the same the states possess very different symmetries as evident by comparing first two rows of figure 3 with the third and fourth row. This clearly shows that the direction of the magnetic field can play a very significant role on the structure of the steady states even though it may not influence their number.

In general, under the presence of the isotropic magnetic field the symmetries subspace merge and the nonequilibrium steady state changes accordingly. On the other hand, the magnitude of the magnetic field does not play a relevant role in the form of the NESS. To make this claim concrete, we compare the NESSs with and without a magnetic field using the fidelity, defined as

$$\begin{equation}\mathcal{F}\left({\rho }_{1},{\rho }_{2}\right)=\mathrm{Tr}\sqrt{\sqrt{{\rho }_{1}}{\rho }_{2}\sqrt{{\rho }_{1}}}.\end{equation} \tag{ 5 }$$

In figure 4 we display the fidelity of the square system with edge coupling to the baths. When the magnitude of the magnetic field is B = nπ, with n integer, the fidelity is maximum as the RDMs are the same. For a finite but not zero value of B there is a discontinuous change in the value of the fidelity due to the symmetry breaking of the Liouvillian. Even after breaking the symmetry the value of the fidelity is always close to one, as the RDMs are similar. This shows that the NESS is robust under the presence of the magnetic field. The main difference appears when B = π/2, but even in this case the fidelity is almost 0.99. In case of the single plaquette face connected system the fidelity always remains at unity, indicating that the NESS is invariant under the application of a magnetic field.

Zoom In Zoom Out Reset image size

For the cube scenario the fidelity shows a slightly more complicated behaviour as shown in figure 5. In the edge-coupling case (panel (a)) the difference in fidelity between having no field and a non-zero field is bigger than in the single plaquette case, but it is still very small ≈1 [similar observations holds true for face coupled case, panel (b)]. Besides, the fidelity presents several local minima, all of them only slightly lower than the flat value. Thus, the NESS is changed because once there is a magnetic field it does not have well defined symmetries, but it is still very similar to the non-magnetic field case implying that the magnitude of the field itself plays a minimal role.

Zoom In Zoom Out Reset image size

For completeness, we have also analyzed a stacked case, i.e. a cubic network with no connections in the z-direction (figure 6). In this case we found 8 different steady states (panels (a)–(p)), where only one is an NESS (panels (a) and (b)), 5 others are pure-states and there are also 2 non-physical RDMs (panels (g) and (h)), meaning that they have zero trace. The existence of these non-physical solutions to a degenerated Lindblad master equation was hypothesised in reference [1] but no examples had been found until now. The bottom two rows correspond the presence of an isotropic or anisotropic field $\to {B}=\left(0,0,B\right)$ or $\to {B}=\left(B,B,B\right)$. Since the vertical links are absent in the stacked case the x- and y-direction magnetic fields play no role, since these only affect the vertical couplings. In this case again we only have one NESS (panels (q) and (r)) whereas all the others are pure states. The pure states show a pattern of symmetries in the simple cubic lattice far more complex than the single plaquette (square) lattice shown in figure 1.

Zoom In Zoom Out Reset image size

To summarize this section, the main idea is to study simple models that are analytically tractable. We chose to present the RDMs visually instead of presenting the complicated analytic solutions obtained using the method described in appendix B. In these simple cases, we evaluated the initial condition independent steady state and found that breaking closed-system symmetries affects the open system properties. Moreover, the RDMs are affected mainly by the direction of the magnetic field and not by its magnitude (see figures 4 and 5: Fidelity deviations from the flat line when a magnetic field is applied). The face coupled system elucidated the role of the system-bath coupling operator and showed that the presence of magnetic field can give rise to non-topological symmetry. Here, even though the presence of magnetic field breaks open-system symmetries it does not affect all RDMs equally. In this case, we observed that the NESS remains invariant even though the other pure-state RDMs were affected. Our simple model also allowed us to find non-physical states which could be identified in the stacked two-plaquette structure with edge coupling to the baths. Thus, our simple model helps us to unravel some of the important consequences of open-system symmetries.

Increasing in the network size greatly increases the complexity making it impossible to obtain analytic results beyond the two plaquette system. Hence, we address this issue numerically by calculating the number of steady states. It is worth noting that even obtaining the initial condition independent steady state density matrices numerically is highly non-trivial task for a degenerated system since the numerical output can be any linear combination which could violate positivity (see appendix B).

Here, we consider three cases of a cubic system illustrated in figure 7. The stacked system consists of N disconnected plaquettes with each plaquette connected to the reservoir via an edge coupling as shown in figure 7(a) (inset). The ladder systems have inter-plaquette connections but a dichotomy exists depending on the system-reservoir coupling. The reservoir can either connect to the faces of the first and last plaquette of the ladder [inset panel (b)], i.e., face connected ladder, or can connect via the edges [inset panel (c)], i.e., edge connected ladder. In the absence of magnetic field the system can have weak and strong symmetries [3] that are analytically intractable. Hence to gain insight into the complex behaviour of the number of symmetry operators on the number of plaquettes in the system we plot the number of steady states (SS) in the top row of figure 7.

Zoom In Zoom Out Reset image size

In the edge connected stacked system, without the constraint of the inter-plaquette connections, the plaquettes are free to permute. This freedom leads to a faster than exponential growth (according to Sterling approximation) in the number of symmetries that is reflected in the fast growing number of steady states as seen in panel (a) (red solid line) for zero magnetic field $\to {B}=\left(0,0,0\right)$. The intra-plaquette symmetry, i.e., flip along the diagonal of a plaquette as discussed previously, plays a rather insignificant role as compared to the symmetry generation due to inter-plaquette permutations. Hence, even when we break the intra-plaquette symmetry by applying an anistropic magnetic field along the z-direction $\to {B}=\left(0,0,B\right)$ [blue dotted line in panel (a)] we still observe an exponential-like growth with slightly smaller number of steady states. Here due to the absence of intra-plaquette connections the isotropic magnetic field scenario is equivalent to the ansiotropic case (blue dotted and green dashed lines overlap). The presence of multiple steady states also affects the nonequilibrium observables including the steady-state particle currents (${J}_{\text{L}}=\mathrm{Tr}\left[\left({A}_{\text{L}1}^{{\dagger}}{A}_{\text{L}1}-{A}_{\text{L}2}^{{\dagger}}{A}_{\text{L}2}\right){\rho }_{\text{ss}}\right]$) [2] as shown in the bottom row of figure 7. Since the direction of the magnetic field affects the symmetries, it also affects the currents in the system with the minimum current obtained when we have maximum number of steady states $\to {B}=\left(0,0,0\right)$. This is simply because the initial condition is a complete dark state of the system ρ(0) ∝ ∑_n=2,4,⋯|n⟩⟨n|−|n⟩⟨n + 2| − |n + 2⟩⟨n| under intra-plaquette flip symmetry which is an eigenstate only when there is no applied magnetic field. In this case only the presence or absence of the magnetic field affects the transport creating a switch that can turn the current on or off.

The face connected ladder system shows a linear dependence between the number of symmetries and the number of plaquettes N. Since the system is placed in a dissipative gradient the inter-plaquette permutations are forbidden and do not contribute to multiple steady states. Hence, in this case the main contribution to the number of steady states comes from the intra-plaquette exchanges. The intra-plaquette symmetry that matters is the 180° rotation, which is preserved when an anisotropic field $\to {B}=\left(0,0,B\right)$ is applied to the system. Hence, as seen in panel (b), the number of steady states remains invariant in the zero field (red solid line) to anisotropic field (blue dotted line) case. When an isotropic magnetic field is applied in this case all symmetries are broken and the system has only one unique steady state (green dashed line). In this case, the current obeys Fouriers law due to the nonlinear coupling with the reservoir even though the system is ballistic [14]. As compared to the edge connected stacked system, where the mere presence of a magnetic field produced nonzero currents, here the direction of the magnetic field strongly affects the transport properties and only an isotropic field is able to produce finite currents.

The edge connected ladder system displays a highly complex behavior with stepwise⁷ increase in the number of steady states as a function of the number of plaquettes figure 7(c). The stepwise behavior is present for an anisotropic field or in the absence of a field, but not in an isotropic field. Hence, the emergence of this behavior seems to be a result of inter-plaquette symmetries. The step-like behavior causes the currents to fluctuate, but as the number of plaquettes N increase the anisotropic and isotropic magnetic field show the same behaviour. Overall, in the three systems in figure 7, different symmetries are broken as a result of applying an anisotropic and isotropic magnetic field. This leads to different manipulation mechanisms of nonequilibrium currents with the presence (edge connected stacked) and direction (face connected ladder) of the magnetic field along with the number of plaquettes (edge connected ladder) playing important controls.

To summarize, this work proposes the control of nonequilibrium properties of a quantum cubic network using a magnetic field. Unlike traditional approaches such as the Aharonov–Bohm effect where the magnitude of the magnetic field is used to control transport [33, 35], we provide a fresh perspective to utilize the direction of the magnetic field. Specifically, this paper shows how a simple tight-binding lattice with engineered couplings to leads could possess symmetries which can be broken systematically by the direction of the applied magnetic field. The symmetry breaking controls not only the steady states but also observables like the currents. Our results do not depend on the choice of parameters such as the magnitude of the magnetic field, coupling strength to the environment, or the specific values of the hopping and on-site potentials as long as the symmetries are maintained. Although this paper focuses on the transport in simple cubic networks, one can envision similar implementations in optics [36] or in new materials such as multi-layer graphene [37] that would extend our work to systems with hexagonal symmetries. We expect our work to trigger significant developments in understanding the connections between various lattice symmetries and open-system symmetries in order to manipulate quantum transport at the nanoscale. A promising research direction would be to utilize the presence of pure states to store energy and strategically release this energy using magnetic field controls to form quantum batteries [38].

This research was supported by the Institute for Basic Science in Korea (IBS-R024-Y2), DM acknowledges the Spanish Ministry and the Agencia Española de Investigación (AEI) for financial support under grant FIS2017-84256-P (FEDER funds), and JC acknowledges support from NSF Grant 1800301 and 1836913.

In this section, we will show how the general many particle representation could be reduced to the single-particle picture. The system Hamiltonain as given by equation (1) reads,

$$\begin{equation}H=-\sum _{r}{\varepsilon }_{r}{c}_{r}^{{\dagger}}{c}_{r}-\sum _{r,d=1,2,3}{t}_{d}{U}^{d}\left(\hat{r}\right){c}_{r+d}^{{\dagger}}{c}_{r}+\mathrm{h}.\mathrm{c}..\end{equation} \tag{ A.1 }$$

The tight-binding system Hamiltonian can easily be expressed in the single particle Fock basis as,

$$\begin{equation}{H}^{1}=-\sum _{r}{\varepsilon }_{r}\vert r\rangle \langle r\vert -\sum _{r,d=1,2,3}{t}_{d}{U}^{d}\left(\hat{r}\right)\vert r+d\rangle \langle r\vert +\mathrm{h}.\mathrm{c}..\end{equation} \tag{ A.2 }$$

Above, the state |r⟩ ≡ |00⋯1_r0⋯⟩ is a single particle Fock state with one particle being present at the rth site. The Lindblad operators connect the single-particle and vacuum state manifolds and are given by,

$$\begin{equation}{A}_{\alpha 1}\equiv {A}_{\alpha 2}^{{\dagger}}=\sum _{n}{c}_{n}^{{\dagger}}{{\Pi}}_{i\ne n}{c}_{i}{c}_{i}^{{\dagger}},\end{equation} \tag{ A.3 }$$

with the sum over n being restricted to n = 1, 5, 9, ⋯ (n = 1, ..., 4) for α = L and n = 3, 7, 11, ⋯ (n = N − 3, ..., N) for α = R for edge (face) connected baths. In the Fock basis representation the Lindblad operator ${A}_{\alpha 1}^{\left(\mathrm{SE}\right)}=\vert n\rangle \langle 0\vert$ with |0⟩ representing the ground state. The form of the Lindbladians does not increase the system Hilbert space exponentially but rather only adds the vacuum state to the single-particle state space. In other words, the dimension of the system is N + 1, N being the number of sites, instead of single-particle dimension N. Thus, we take our invariant basis set (N + 1 states) comprising of the single-particle basis and the vacuum state to construct a single excitation Hamiltonian H^(SE) = H¹ + ɛ₀|0⟩⟨0| with ɛ₀ being the vacuum state energy. The Hamiltonian ensures that there is at most one excitation/particle in the system. In realistic systems, this could be a result of being at low temperatures wherein the probability of finding more than one excitation is very low. Therefore, the quantum master equation governing the single excitation reduced density matrix reads

$$\begin{equation}\frac{\mathrm{d}{\rho }^{\left(\mathrm{S}\mathrm{E}\right)}}{\mathrm{d}t}=-\mathrm{i}\left[{H}^{\left(\mathrm{S}\mathrm{E}\right)},{\rho }^{\left(\mathrm{S}\mathrm{E}\right)}\right]+\sum _{\genfrac{}{}{0pt}{}{k=1,2}{\alpha =\text{L},\text{R}}}{{\Gamma}}_{\alpha k}\mathcal{D}\left[{A}_{\alpha k}^{\left(\mathrm{S}\mathrm{E}\right)}\right]{\rho }^{\left(\mathrm{S}\mathrm{E}\right)}.\end{equation} \tag{ A.4 }$$

Since we are interested in single particle currents and symmetries of the density matrices, dealing with the single excitation reduced density matrix reduces the computation drastically allowing us to deal with large number of sites.

In the simple cubic structures we are tacking we have a quantum Liouvillian $\mathcal{L}$ with several zero eigenvalues. These eigenvalues are a consequence of the existence of one or several non-trivial symmetry operators, π_i, such that they commute with all the generators of the dynamics in the form [π_i, H] = [π_i, A_α,k] = 0  ∀ (α, k) (strong symmetry) or because of the existence of one or several superoperators, Π_i, that commutes with the full Liouvillian operator (weak symmetry) [1]. Finding the steady states of such Liouvillians is highly non-trivial, as by direct diagonalisation we can obtain eigenvectors of the Liouvillian that do not fulfill the requirement for being density matrices (semi-definite positivity). Only if we know a complete set of operators, i.e., the set {π_i, H, A_α,k} of all jump and symmetry operators that generate the entire algebra of the operators Hilbert space, we can calculate all the density matrices by block diagonalising the Liouvillian matrix.

In the worst-case scenario we have a Liouvillian with N zero eigenvalues, but no information about the symmetries. This Liouvillian could be reduced to a smaller one if some but not all of the symmetries are known. As the degeneracy comes from the symmetry we know that there are N Hermitian matrices ρ_i, s.t. $\mathcal{L}{\rho }_{i}=0$ (which we will refer as zero-eigenmatrices), belonging to different symmetry subspaces. It is easy to prove that these matrices should be orthogonal in the sense that Tr  ρ_i ρ_j = 0 if i ≠ j. Most likely, these matrices would have non-zero trace and they may be semi-definite positive, but there is also the possibility of having matrices with zero trace if they belong to a subspace with mixed symmetries (see references [1, 3] for discussion). We will not consider the latter case below, as it is untypical, but the method proposed here can be trivially extended to a case with zero-trace zero-eigenmatrices. By direct diagonalisation of the Liouvillian matrix we can obtain N eigenmatrices, ${\tilde {\rho }}_{i}$, s.t. $\mathcal{L}{\tilde {\rho }}_{i}=0$. These matrices need not be Hermitian, semi-definite positive or have unit trace as they can be just a combination of the real density matrices, i.e., they do not belong to a specific symmetry subspace. Furthermore, they do not necessarily fulfill the orthogonality condition. The problem thus can be summarized as follow: we have a set of N zero-eigenmatrices $\left\{{\tilde {\rho }}_{i}\right\}$ that are a linear combination of N density matrices ρ_i, how can we recover the set {ρ_i} from $\left\{{\tilde {\rho }}_{i}\right\}$?

First, we need the matrices to be Hermitian. This can be done by generating new matrices ${\rho }_{i}^{\text{H}}={\tilde {\rho }}_{i}+{\tilde {\rho }}_{i}^{{\dagger}}$. Besides, we need the matrices to form an orthogonal set. We then apply a Gramm–Schmidt type algorithm in the following form;

$$\begin{align}\hfill {\rho }_{1}^{\text{o}}& ={\rho }_{1}^{\text{H}}.\hfill \\ \hfill {\rho }_{2}^{\text{o}}& ={\rho }_{2}^{\text{H}}-\frac{\mathrm{Tr}{\rho }_{1}^{\text{o}}\enspace {\rho }_{2}^{\text{H}}}{\mathrm{Tr}{\rho }_{1}^{\text{o}}\enspace {\rho }_{1}^{\text{o}}}{\rho }_{1}^{\text{o}}\hfill \\ \hfill {\rho }_{3}^{\text{o}}& ={\rho }_{3}^{\text{H}}-\frac{\mathrm{Tr}{\rho }_{1}^{\text{o}}\enspace {\rho }_{3}^{\text{H}}}{\mathrm{Tr}{\rho }_{1}^{\text{o}}\enspace {\rho }_{1}^{\text{o}}}{\rho }_{1}^{\text{o}}-\frac{\mathrm{Tr}{\rho }_{2}^{\text{o}}\enspace {\rho }_{3}^{\text{H}}}{\mathrm{Tr}{\rho }_{2}^{\text{o}}\enspace {\rho }_{2}^{\text{o}}}{\rho }_{2}^{\text{o}}\hfill \\ \hfill & {\vdots}\hfill \\ \hfill {\rho }_{N}^{\text{o}}& ={\rho }_{N}^{H}-\sum _{j=1}^{N-1}\frac{\mathrm{Tr}{\rho }_{j}^{\text{o}}\enspace {\rho }_{N}^{\text{H}}}{\mathrm{Tr}{\rho }_{j}^{\text{o}}\enspace {\rho }_{j}^{\text{o}}}{\rho }_{j}^{\text{o}}.\hfill \end{align} \tag{ B.1 }$$

The set $\left\{{\rho }_{i}^{\text{o}}\right\}$ is now a set of orthogonal zero-eigenmatrices of the Liouvillian. The only problem remaining is that this set may contain non-semi-definite positive matrices, meaning that the eigenvalues of some of the elements of $\left\{{\rho }_{i}^{\text{o}}\right\}$ may be negative. To transform them in a set of positive zero-eigenmatrices and preserving the orthogonality of the system we may apply a rotation in the N real space. We define a column vector of the zero-eigenmatrices ${\to {\rho }}^{\enspace 0}{:=}{\left({\rho }_{1}^{\text{o}},\enspace {\rho }_{2}^{\text{o}},\dots ,{\rho }_{N}^{\text{o}}\right)}^{\text{T}}$ and the rotation is given by a unitary matrix U_N(θ₁, θ₂, ..., θ_k) with {θ₁, θ₂, ..., θ_k} being a set of Euler angles with $k=\frac{{N}^{2}-N}{2}$. For a fixed set of angles we can obtain a new set of rotated eigenmatrices as $\to {\rho }\left({\theta }_{1},{\theta }_{2},\dots ,{\theta }_{k}\right)={U}_{N}\left({\theta }_{1},{\theta }_{2},\dots ,{\theta }_{k}\right){\to {\rho }}^{\enspace 0}$.

To fix the parameters to find the rotation that makes all the density matrices positive we define a function

$$\begin{equation}F\left({\theta }_{1},{\theta }_{2},\dots ,{\theta }_{k}\right)=\sum _{i=1}^{N}\sum _{j=1}^{d}{\lambda }_{i,j}\left({\theta }_{1},{\theta }_{2},\dots ,{\theta }_{k}\right)-\left\vert {\lambda }_{i,j}\left({\theta }_{1},{\theta }_{2},\dots ,{\theta }_{k}\right)\right\vert ,\end{equation} \tag{ B.2 }$$

with d being the dimension of the density matrices space and λ_i,j(θ₁, θ₂, ..., θ_k) being the jth eigenvalue of the ith matrix from the vector $\to {\rho }\left({\theta }_{1},{\theta }_{2},\dots ,{\theta }_{k}\right)$. It is clear that when the vector $\to {\rho }\left({\theta }_{1},{\theta }_{2},\dots ,{\theta }_{k}\right)$ contains only positive matrices if the corresponding function F(θ₁, θ₂, ..., θ_k) is equal to zero. Therefore, we can solve the equation $F\left({\theta }_{1}^{{\ast}},{\theta }_{2}^{{\ast}},\dots ,{\theta }_{k}^{{\ast}}\right)=0$, or maximize the function F, in order to find the angles that transform our set into a set of positive matrices.

$$\begin{equation}{\to {\rho }}^{\text{P}}=U\left({\theta }_{1}^{{\ast}},{\theta }_{2}^{{\ast}},\dots ,{\theta }_{k}^{{\ast}}\right){\to {\rho }}^{\enspace 0}.\end{equation} \tag{ B.3 }$$

Finally, we normalize each matrix in the vector ${\to {\rho }}^{\text{P}}$ by doing ${\rho }_{i}={\rho }_{i}^{\text{P}}/\mathrm{Tr}{\rho }_{i}^{\text{P}}$ obtaining the desired set of density matrices.

By the use of this method we can calculate the steady-states of any degenerated Liouvillian without knowing the symmetry operators. The main difficulty of this method is to perform a minimization over $\frac{{N}^{2}-N}{2}$ angles, being N the degeneracy number, and the calculation of the function F that requires the diagonalisation of N matrices. For big, highly degenerated systems the computational problem is very hard and the use of efficient numerical algorithms for minimisation and diagonalisation is required. For small systems with a low degeneration it is possible to find the steady-states, as we did with all the plots of this paper. Finally, if some but not all of the symmetry operators are known it is possible to reduce the complexity of the problem by applying the orthonormalisation algorithm only to the symmetry eigenspaces.