Magnetic field induced symmetry breaking in nonequilibrium quantum networks

We study the effect of magnetic field on the nonequilibrium transport properties of a general cubic quantum network described by a tight-binding Hamiltonian. We demonstrate that the symmetry of open systems can be manipulated by the direction of the magnetic field. Starting with all the symmetries preserved in absence of a field, the anisotropic and isotropic field systematically break the symmetries, influencing all nonequilibrium properties. For square and cubic structure, we are able to identify the steady states that comprise of pure states, bath-dependent states (nonequilibrium steady states), and also nonphysical states. As an application, we show numerically for large cubic networks that the symmetry breaking can control nonequilibrium currents and that different environmental interactions can lead to novel features which can be engineered in artificial super-lattices and cold atoms.

Introduction.-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]. 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.
Many theoretical consequences of symmetry have been studied [6], but little is known about the connections between structural symmetries and open quantum system symmetries. 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 systems with structural symmetries found in many natural materials [7]. Recently, several experimental setups have been made accesible 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][11][12].
On the theoretical side, transport in quantum cubic networks have been studied using different approaches. In one-dimensional linear systems made of spins [13,14] and harmonic oscillators [15] the ballistic transport is suppressed when a noise channel is coupled to the system. In a square spin lattice of higher dimensions there is a richer behaviour and both ballistic and diffusive subspaces can coexist [16,17]. In higher dimensions it has been proved that quantum walks in hypercubes hit exponentially faster in comparison to their classical counterparts even if these hypercubes are distorted [18][19][20].
Inspired by these results and the structural simplicity of cubic networks in this letter we aim to address if the underlying symmetry can help us obtain the steady states. We analyze the square and cubic networks driven far from equilibrium and obtain the underlying open system symmetry operators and consequently the steady state system density matrices that reflect the symmetries. Furthermore, we study cubic networks and show that using Abelian magnetic fluxes [21] the symmetries in these networks can be systematically broken and can be harvested to control nonequilibrium transport properties. In particular, we find that the direction of the magnetic field plays a crucial role to control transport properties. Cubic Network in a Magnetic Field.-The system of interest is a slab of 3D cubic network of non-interacting fermions subjected to an external magnetic field illustrated in Fig. 2. The tight-binding representation of the system Hamiltonian is given by, 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 [22] U 1 (r) = 1, U 2 (r) = exp i2πB x − y − 1 2 , and U 3 (r) = exp [−iπB (x − y)] with the cube length a = 1. The phases are obtained under the Hasegawa gauge [23,24] that leads to a translationally invariant Hamiltonian [25][26][27].
The above magnetic field descriptions have been experimentally realized in artificial super-lattices [28][29][30][31] and cold atoms [21,32,33]. The latter have a high level of control on the parameters and hence form the ideal test-bed. Our system is connected to two reservoirs as shown in Fig. 2 allowing us to explore the effects of magnetic fields on spatially symmetric nonequilibrium dissipative systems. We engineer the system-bath coupling 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 [34,35]. The fermionic operators d, d † correspond to the annihilation and creation operators of the reservoir H α = k k d † k,α d k,α with ν k being the coupling strength of each of the free fermions of the reservoir to the system. The specific form of the coupling allows the reservoirs to locally inject or extract a particle between the single particle states corresponding to the sites connected to the reservoir and the ground state. The dynamics of the system reduced density matrix is governed by a Markovian dissipative Lindblad quantum master equation [36][37][38][39] that reads, Here ) allow local injection (k = 1) or extraction (k = 2) of particles into the system H [Eq. 3]. 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 injec- being the Fermi-Dirac distribution. The ratio of the rates Γ α1 /Γ α2 = exp[−β α 0 ] obeys local-detailed balance where β α = 1/k B T α represents the inverse temperature of the α-th bath and 0 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,2,4,6]. Square and Cube.-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 even obtain all the the  [40] for a single plaquette (square) system. Panels a and b illustrate the edge coupled plaquette (top view) with sites labeled (1-4) and face coupled plaquette (side view) respectively. Next two rows correspond to the edge coupled plaquette (panels c-h) without (panels c-d and g-h) and with (panels e-f) an anisotropic magnetic field in the z direction. Panels i-t correspond to the face coupled plaquette with (panels i-j, m-n, and q-r) and without (panels k-l, o-p, and s-t) an anisotropic magnetic field. The first and the third columns correspond to the real parts of the density matrices, whereas second and fourth columns are the imaginary parts. Each panel is normalized by the largest value and is the density matrix plotted in the single particle basis i (i = 1, · · · , 4) with 0 being the vacuum state of the system. The parameters used for the calculations are: ε1 = ε3 = 1, ε2 = ε4 = 1.  Figure 1 shows the various steady states for a single plaquette (square) system illustrated in panels a and b. The edge coupled (panel a) and face coupled (panel b) reservoir posess different open system symmetries, giving rise to different steady states. In case of edge coupled plaquette, the open system symmetry that survives in absence of magnetic field is the flip symmetry along the diagonal (indicated by the dashed line in Fig. 1a). The corresponding unitary symmetry operator is U = exp (|2 4| + |4 2|). The symmetry operator block diagonalizes the Liouvillian giving rise to two steady states as indicated in panels c-d and g-h, 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 c-d) and the other is a pure state independent of the reservoir properties and depending only on the system symmetries (panels g-h). The pure state turns out to be anti-symmetric under the exchange of states |2 and |4 as seen in panel g. Such states carry zero current [42] and are also referred to as dark or a compact localized states [43]. In presence of an anisotropic (z-direction) or isotropic magnetic field (panels e-f) the edge coupled system displays a single NESS that differs from the one obtained in the absence of magnetic field (panels c-d).
The face coupled plaquette in the absence of magnetic field (panels i-j, m-n, and q-r) gives three physical steady states with two pure states that are anti-symmetric with respect to exchange of sites |2 −|4 (panel q) and |1 −|3 (panel m). The corresponding symmetry operators are 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 180 • rotation symmetry still persists. Additionally another symmetry emerges that gives us three steady state solutions even in the presence of a magnetic field (panels k-l, o-p, and s-t). In this case all the states are physical and we still find only one NESS. The two pure states (panels o-p and s-t) 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 system provides a range of unexpected results, like the presence of extra symmetries with the magnetic field, invariance or change in the steady state w.r.t. magnetic field, etc., that help describe the complete nonequilibrium picture. We have also obtained the density matrices for the cubic system (both edge and face coupled), see supplementary material [41]. Increasing a single plaquette, i.e., cubic network, increases the complexity drastically. The edge coupled stacked cube (cube without links between the plaquettes in the z direction) possess 6 steady states for B = (0, 0, 0), 5 steady states for B = (0, 0, B) (anisotropic field), and 1 for B = (B, B, B) (isotropic field). In all cases only one steady state is a genuine NESS in the sense that it depends on the bath parameters. All the other states are pure states and are eigenstates of the system Hamiltonian. 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. This is discussed in the supplementary material by calculating the fidelity between the NESS with B = (0, 0, 0) and the NESS with B = (B, B, B) [41]. For completeness, we have also analised a stacked case, i.e. a cubic network with no connections in the z-direction. In this case we found 8 different steady states, where only one being a NESS, 5 others are pure-states and there are also 2 non-physical density matrices, meaning that they have zero trace. The existence of these non-physical solutions to a degenerated Lindblad master equation was hypothesised in Ref. [1] but no examples had been found until now. Cubic Slab.-Increase in system size leads to a highly complex symmetry problem despite the simple geometrical network. Beyond the two plaquette system it becomes impossible to analytically obtain the steady states and hence we address the issue by counting the number of steady states numerically. Obtaining the steady state density matrices, is highly non-trivial numerically, since the numerical output can be any linear combination which could violate positivity (See supplementary material [41] Sec. II). But, the number of steady states, including both physical and non physical, can be easily obtained by counting the number of zero eigenvalues of the dissipative generator.
Here, we consider three cases of a cubic system illustrated in Fig. 2. The stacked system consists of N disconnected plaquettes with each plaquette connected to the reservoir via an edge coupling as shown in the inset of panel a. The ladder systems have inter-plaquette connections but a dichotomy exists depending on the systemreservoir 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 [2] 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 in the top row of Fig. 2 [44].
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 B = (0, 0, 0). 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 B = (0, 0, B) [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 cur- [2] as shown in the bottom row of Fig. 2. 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 B = (0, 0, 0). This is simply because the initial condition is a complete dark state of the system ρ(0) = 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 . Intuitively, 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 B = (0, 0, B) 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 [17]. 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 [45] increase in the num-ber of steady states as a function of the number of plaquettes Fig. 2c. The stepwise behavior that is present without a magnetic field and with an anisotropic field. Hence, the emergence of this behavior seems to be a result of inter-plaquette symmetries, although intra-plauette symmetries play a crucial role in this system causing the number of steady states to differ from the zero magnetic field case. The step-like behavior also causes the currents to fluctuate, but as the number of plaquettes N increase the currents merge in the presence of a magnetic field. Overall, in the three systems in Fig. 2, 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.
Conclusions.-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 [42,46], we provide a fresh perspective to utilize the direction of the magnetic field. We specified how a trivial 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. Although we dealt with simple cubic networks, one can envision similar implementations in optics [47] using an array of nanoparticles or in new materials such as multi-layer graphene [48] 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 [49].
In the Fock basis representation the Lindblad operator A (SE) α1 = |n 0| with |0 representing the ground state. The form of the Lindbladians doesn't 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 1 + ε 0 |0 0| with ε 0 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 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.

CALCULATION OF THE DENSITY MATRICES OF A DEGENERATED LIOUVILLIAN
In the simple cubic structures we are tacking we have a quantum Liouvillian 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 symmetry operators, in the sense that the set {π i , H, A α,k } including all the jump and symmetry operators that generate the entire algebra of operators of the system, 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 Liovillian could be the 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. Lρ 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 Refs [1,2] for discussion). We will not consider the latter case below, as it is an untypical case, 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,ρ i , s.t. Lρ 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 summarised as follow: We have a set of N zero-eigenmatrices {ρ i } that are a linear combination of N density matrices ρ i , how can we recover the set {ρ i } from {ρ i }? First, we need the matrices to be Hermitian. This can be done by generating new matrices ρ H i =ρ i +ρ † i . Besides, we need the matrices to form an orthogonal set. We then apply a Gramm-Schmidt type algorithm in the following form; The set {ρ o i } 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 {ρ o i } 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 ρ 0 := (ρ o 1 , ρ o 2 , · · · , ρ o N ) T and the rotation is given by a unitary matrix U N (θ 1 , θ 2 , · · · , θ k ) with {θ 1 , θ 2 , · · · , θ k } being a set of Euler angles with k = N 2 −N
To fix the parameters to find the rotation that makes all the density matrices positive we define a function with d being the dimension of the density matrices space and λ i,j (θ 1 , θ 2 , · · · , θ k ) being the jth eigenvalue of the ith matrix from the vector ρ(θ 1 , θ 2 , · · · , θ k ). It is clear that when the vector ρ(θ 1 , θ 2 , · · · , θ k ) contains only positive matrices if the corresponding function F (θ 1 , θ 2 , · · · , θ k ) is equal to zero. Therefore, we can solve the equation F (θ * 1 , θ * 2 , · · · , θ * k ) = 0, or maximize the function F , in order to find the angles that transform our set into a set of positive matrices.
Finally, we normalise each matrix in the vector ρ P by doing ρ i = ρ P i /Tr ρ P i obtaining the desired set of density matrices.

EFFECT OF THE MAGNETIC FIELD ON THE NESS
As discussed in the main text, the main effect of the magnetic field is to break the degeneracy of the Liouvillian, giving rise to a reduced number of steady states. This happens for any finite value of the magnetic field. In this section, we discuss the effect of the magnitude of the magnetic field on the form of the density matrices. To compare different density matrices we have used the Fidelity, defined as In Fig. 3 we display the Fidelity between the NonEquilibrium Steady States (NESSs) of the 1-plaquette system with and without magnetic field 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 density matrices 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 very close to one, as the density matrices are very similar. This shows that the NESS is robust under the presence of the magnetic field. The main difference appear when B = π/2, but even in this case the fidelity is almost 0.99. In case of the 1-plaquette face connected system the fidelity always remains at unity, indicating that the NESS is invariant under the application of a magnetic field. For the 2-plaquette scenario the fidelity shows a slightly more complicated behaviour as it is shown in Figures 4 (edge coupling) and 5 (face coupling). In the edgecoupling case the difference between having no field and a very small one is bigger than in the 1-plaquette case, but it is still very small. Besides, the fidelity presents several local minima, all of them only slightly lower than the flat value.
We can conclude from this analysis that the main effect of the magnetic field acting on the system is to break the degeneracies. 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. Thus, implying that the magnitude of the field itself plays a very minimal role.

STEADY STATES OF A CUBE
We study the effects of magnetic field on a system comprised of two plaquettes (cube) connected in different ways to the reservoir. Figure 6 looks at the stacked system, i.e., cube without vertical links, connected via the edges to the two reservoirs. The first two rows of Fig. 6 are the density matrices in the absence of a magnetic field. Using the procedure outlined in Sec. we evaluated all the density matrices. Surprisingly, we find that one of the density matrices is unphysical (panels gh) since it has trace zero and it contains both positive and negative populations in the Fock basis. The existence of such density matrices was proposed in Ref. 1, but till date no concrete examples were found. Moreover, as explained in the main text all the density matrices are pure states except one (panels a-b) which is a NESS that depends on the properties of the two reservoirs. The bottom two rows correspond the presence of an isotropic or anisotropic field B = (0, 0, B) or B = (B, B, B). Since the vertical links are absent in the stacked case the xand 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-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 illustrated in the main text.
Next we look at the face connected cube in Fig. 7. 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 B = (0, 0, B) (panels k-t) and only 1 state for an isotropic field B = (B, B, B) (panels u-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 Fig. 7 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 drastically. Unlike the single plaquette case described in the main text the NESS (panels a-b, k-l, and u-v) is not invariant when the magnetic field is turned on.
Lastly, we look at the cubic lattice connected on the edge with reservoirs in Fig. 8. Without magnetic field we obtain 6 steady states, which reduce to 5 when an anisotropic field is applied, which further reduces to 1 in presence of an isotropic field. Similar to all previous cases, we have only 1 NESS for each direction of the magnetic field and the rest are all pure states. The presence of an isotropic field destroys all open system symmetries and gives a unique steady state. In all cases for the two plaquette system, all the pure states in the absence of the magnetic field are always real which obtain imaginary parts as soon as a magnetic field is applied. As evident from all the cases above the symmetries of the reduced density matrix are non trivial and exhibit complex behaviour even for a simple two plaquette structure. * jythingna@ibs.re.kr † manzano@onsager.ugr.es ‡ jianshu@mit.edu [1] B. Buča and T. Prosen, New J. Phys. 14, 073007 (2012).