Ultrafast spin dynamics in inhomogeneous systems: a density-matrix approach applied to Co/Cu interfaces

Ultrafast spin dynamics on femto- to picosecond timescales is simulated within a density-operator approach for a Co/Cu bilayer. The electronic structure is represented in a tight-binding form; during the evolution of the density operator, optical excitation by a femtosecond laser pulse, coupling to a bosonic bath as well as dephasing are taken into account. Our simulations corroborate the importance of interfaces for ultrafast transport phenomena and demagnetisation processes. Moreover, we establish a reflow from Cu $d$ orbitals across the interface into Co $d$ orbitals, which shows up prominently in the mean occupation numbers. On top of this, this refilling manifests itself as a minority-spin current proceeding several layers into the Cu region. The present study suggests that the approach captures essential ultrafast phenomena and provides insight into microscopic processes.


Introduction
The field of spin-electronics comprises a variety of effects that often show up in inhomogeneous systems and as responses to a static cause. In recent years, however, the attention has shifted to ultrafast processes which, for example, are pushed by femtosecond laser pulses or by electromagnetic terahertz radiation. Understanding the fundamentals of ultrafast spin dynamics requires theories which allow to identify the essential mechanisms for, e. g., switching magnetisation or generating spin-polarised currents, provide intuitive but detailed pictures of these processes, and are able to reproduce experimental results.
Ultrafast spin dynamics in inhomogeneous systems covers time scales from femtoto picoseconds. On the femtosecond timescale time-dependent density-functional theory (TDDFT) offers the perhaps most fundamental and detailed description [1]. Since such computations are very demanding, they are presently limited to a short duration (a few ten femtoseconds) and a small sample size (a few unit cells). On top of this, the electron system is not coupled to thermodynamic baths (e. g., phonons or magnons).
Töws and Pastor developed a many-particle theory of ultrafast demagnetisation in ferromagnetic transition metals [2,3], in which spin-orbit coupling and a laser pulse are considered. On top of this, the intra-atomic Coulomb interaction is taken into account, making analytical or numerical solutions very demanding. An application to Ni showed that demagnetisation effects are short-ranged [3], allowing to use small clusters in configuration space.
On the picosecond timescale, at which thermalisation takes place, the subsystems of electrons, phonons, and magnons may be described by their individual temperatures that are obtained from rate equations, as in the three-temperature model [4]. The detailed geometry of a sample, for example inhomogeneities, is usually regarded as side issue.
The above calls for a versatile theoretical single-particle approach that covers all timescales on equal footing but does not suffer from restrictions on a sample's geometry and size. On top of this, it should account for both the excitation of the electron system by a laser pulse and the coupling of the electron system to thermodynamic baths in order to achieve thermalisation. In this paper we report on a theoretical method which aims at fulfilling the previously mentioned objectives. The starting basis is a tight-binding description of the electronic structure in configuration space which is not only versatile but as well provides an intuitive interpretation of the results. As in TDDFT simulations, we obtain mean occupation numbers and currents by evolving the one-particle density operator in time. The occupation numbers and currents are fully resolved with respect to time, site, orbital, and spin, thereby allowing a detailed analysis. Thermalisation is achieved by coupling the electron system to a bosonic bath. As an application we choose spin dynamics at Co/Cu interfaces which due to their epitaxial growth serve as paradigm system in spin-electronics.
Some theories of transport in a ferromagnetic/nonmagnetic sample concentrate on majority-spin electrons, for example superdiffusive spin transport [5,6]. Strong asymmetries of spin-dependent relaxation times and transition probabilities explain that after the laser excitation mostly hot majority-spin electrons traverse the interface and propagate rapidly from the ferromagnetic into the nonmagnetic region of the sample. This scenario is corroborated by experimental results on Co/Cu systems that show that demagnetisation is initially governed by spin-polarized transport starting at the interface [7] (for Fe/Au see Refs. [8] and [9]). The above sequence of events is complemented by a backflow mechanism [10,11]: empty minority-spin Co d orbitals become occupied by Cu d electrons from sites adjacent to the Co/Cu interface. The Co d orbitals can be located above or below the chemical potential. In the latter case they have to be depleted by a laser excitation before being refilled. Chen et al explained this 'reflow' by resonant optical excitation by the 1.5 eV laser pump pulse; it is in line with the high density of states of both Co and Cu in the relevant energy range. The reflow mechanism is reminiscent to optically induced inter-site spin transfer (OISTR; Ref. [12]).
The reflow scenario in Co/Cu is fully confirmed by our simulations. As will be demonstrated in this paper, we identify which orbitals are essential for the above mechanism. Furthermore, the currents triggered by the reflow are investigated for times considerably after the pulse, thereby going beyond the TDDFT simulations by Chen et al [10]. Space-time maps of the associated spin-resolved currents reveal how deep the reflow penetrates into the Cu region of the samples. This paper, whose main purpose is to introduce to the approach and to demonstrate its key capabilities, is organised as follows. In Section 2 we sketch our theoretical approach by addressing sample geometry (Section 2.1), electronic structure (Section 2.2), and evolution of the density operator (Section 2.3). In Section 3 we report on selected features of spin dynamics at Co/Cu interfaces (Section 3.1): occupation numbers and demagnetisation (Section 3.2) as well as spin transfer across the interface (Section 3.3). An outlook is given in Section 4, while further details are comprised in an appendix.

Sample geometry
A sample is represented as a cluster of atoms in configuration space. In order to be flexible, a cluster is built from blocks, each of which contains identical unit cells; the unit cells of the individual blocks may differ from each other. This approach allows to construct samples in one, two or three spatial dimensions, for example chains, heterostructures or homogeneous systems with defects (e. g., vacancies). A typical cluster's shape is a parallelepiped.

Electronic structure
The electronic structure of a beforehand defined cluster is described within the Slater-Koster formulation of the tight-binding (TB) method [13]. For each site (atom) we consider exchange-split s, p, and d orbitals and take into account spin-orbit coupling. The spin quantisation axis is the z axis in the global frame. The tightbinding parameters are either taken from literature or computed from first principles. Furthermore, the boundary conditions (open or closed) are specified for each edge direction of the cluster parallelepiped.
The result of the above setup is the time-independent Hamiltonian matrix H 0 whose eigenvectors enter the evolution of the one-particle density matrix P(t); see Section 2.3 below. Recall that the TB parameters do not depend on time, which is a reasonable assumption as long as any time-dependent perturbation is weak. This implies that a so-called transient band structure [14] can show up in a simulation only as variation of the mean occupation numbers of the eigenstates ofĤ 0 ; the eigenstates and -energies themselves are not affected.
The spin dynamics is analyzed in terms of theĤ 0 eigenstates |n (termed 'eigenstate basis' in this paper) or in the 'site-orbital basis' of site-dependent orbitals |k, α (k site index, α multi-index comprising orbital and spin quantum numbers), as the occasion may require.

Evolution of the density operator
The evolution of the one-particle density operatorρ(t) of the electron system is determined by a laser pulse and by coupling to a bosonic bath.
A laser pulse is specified by its electric-field vector (direction of incidence with respect to the global frame, amplitude, and polarisation), the photon energy, and the shape of a Gaussian envelope function, the latter typically of a few fs width. The pulse induces optical dipole transitions between the orbitals, whose matrix elements are condensed into the perturbation matrix V(t); cf. Ref. [3]. The spin-conserving transitions obey the angular-momentum selection rules δl = ±1 and δm = 0, ±1, depending on the laser's polarisation.
The electron system is thermalised by coupling it to a fictitious bosonic bath, the latter mimicking phonons or magnons. The coupling is described by Lindblad superoperators [15,16] that determine how rapidly the laser-excited electron system will relax toward thermal equilibrium. The Lindblad superoperator of an operatorF acts on the density operatorρ and accounts for energy transfer from the electron system to the bath and vice versa; it describes also dephasing. The number of electrons is conserved since L[F ](ρ) is traceless. A jump operatorF nm ≡ |n m| mediates an inelastic transition from theĤ 0 eigenstate |m with energy ε m into the eigenstate |n with energy ε n , thereby transferring an energy of ∆ε from the electron system into the bath (ε m > ε n ) or vice versa (ε n > ε m ). The Lindbladian L[F nm ](ρ(t)) forF =F nm at time t, follows from Equation (1); ρ nm (t) ≡ n|ρ(t)|m . The total Lindbladian is then obtained by a weighted sum over allL nm . If n = m, the strengths γ nm account for the Bose-Einstein distribution function f BE of the bath as well as for spin-conserving (sc) and spin-flip (sf) transitions. For the latter we introduce transition rates γ sc and γ sf , respectively.
If n = m, the Lindblad operatorsL nn are energy-conserving and do not represent transitions. Instead they cause pure dephasing by reducing the off-diagonal elements of ρ; their diagonal elements vanish. The dephasing rate γ dp determines the timescale on which modulations of the mean occupation numbers are quenched.
Representing the operators as matrices, the evolution of the density matrix P(t) (with matrix elements ρ nm (t)) is given by the Lindblad equation (in Hartree atomic units). P(t) is not diagonal because of the perturbation V(t).
The outcome of the Lindblad equation are mean occupation numbers ρ nn (t) of the eigenstates |n and -after transformation into the site-orbital basis -of the orbitals |k, α , ρ kα,kα (t). The effect of the perturbation and the coupling to the bath is identified as the difference ∆n(t) ≡ ρ nn (t) − ρ (0) nn , that is the deviation of the nonequilibrium occupancy ρ nn (t) with respect to the one in equilibrium (ρ (0) nn , an element of P 0 ). Following a derivation by Mahan [21], the current flowing from site l to site k at time t can be expressed as in which are blocks of the density matrix P(t) and the Hamilton matrix H 0 , respectively, both expressed in the site-orbital basis. The submatrices p σσ lk (t) and h σσ kl are indexed by orbital. For details see Appendix C.
The spin-polarised currents resolved with respect to the µ-th component are given by in which Σ µ is a block Pauli matrix and [·, ·] + is the +-commutator. The z-resolved spin-polarised current studied in this paper reads that is the usual difference of spin-up and spin-down currents. Equation (6) may be construed broadly as follows. The laser pulse produces local imbalances of occupation (that is electric polarisation), which themselves cause siteoff-diagonal ρ terms. The larger the imbalance and the stronger the hopping strength between the sites, the larger the current. Recall that the hopping strengths h kα,lβ determine the bandwidths and, thus, the group velocities of the electrons.

Setup for Co/Cu systems
To illustrate the approach and its capabilities we deliberately choose rather simple Co/Cu systems which have the advantage that numerical results are easily visualised and interpreted but nevertheless exhibit the relevant physical processes.
A two-dimensional cluster of 20 Co and 20 Cu atoms simulates dynamics in a thin fcc Co layer on a Cu substrate with a (001) interface. Periodicity is applied in the directions parallel to the interface (y and z). For studying spin transfer we replaced this cluster by one with 10 Co and 30 Cu atoms.
The tight-binding parameters were taken from Ref. [22] but with the on-site energies of bulk Co and Cu adjusted to achieve aligned Fermi energies. The Co-Cu hopping parameters (across the interface) are taken as arithmetic means of the respective ones of Co and Cu; we are aware that this may be regarded as too rough an approximation but recall that similar heuristic settings worked for tunnel junctions [23].
The electron system is excited by a linearly polarized laser pulse with a photon energy of 1.55 eV and a width of the Gaussian envelope of 10 fs. The Lindblad parameters γ nm defined in Equation (4) are chosen such as to reproduce relaxation and dephasing times in experiments: γ sc = 2 · 10 −4 fs −1 , γ sf = 2 · 10 −6 fs −1 , and γ dp = 5 · 10 −2 fs −1 , which correspond to lengths of times of 5 ps, 500 ps and 20 fs, respectively. The results presented in this paper do not depend crucially on their actual values. The bath temperature is set to T = 300 K in all simulations.

Mean occupation numbers and demagnetisation
As addressed in the introduction (Section 1), temperature models are often used for describing spin dynamics. For example, time-resolved spectroscopic data are approximated by Fermi-Dirac distributions, yielding the evolution of both electron temperature and chemical potential. Such an approach seems reasonable at times after a fs laser pulse, but needs to be verified for times during or short after a laser pulse.
• Long before the laser pulse, the initial thermalisation produces a Fermi-Dirac distribution of the occupation numbers [panel (a) in Figure 1]. ‡ Next, electron states close to the chemical potential become depopulated or populated, but the distribution is changed slightly (not shown); this is explained by the rather weak amplitude of the laser pulse at these times.
• With increasing amplitude, also electron states well below the chemical potential µ become strongly depopulated; simultaneously, hot electrons appear at energies far above µ [panel (b) shows the distribution at the maximum amplitude]. The distribution is modulated in line with the frequency of the pulse (recall that dipole transitions induce both excitations and de-excitations). ‡ Figure 1 comprises selected snapshots of the evolution animated in Figure B1. The laser pulse produces mainly a depopulation of minority-spin electrons below µ and a population of hot electrons. The latter is larger, the larger the maximum amplitude of the pulse. This finding is confirmed qualitatively by occupation-number profiles computed within TDDFT for similar Co/Cu systems.
The pronounced structure within the distribution of minority electrons is analysed further by projecting the eigenstates onto the involved orbitals (s, p, and d; Figure 2 panels a -c). Panel c shows that this structure is composed mainly of eigenstates of d orbital character. While those states with reduced occupation below the chemical potential µ are identified as Cu states, eigenstates with enhanced occupation exactly at and above µ exhibit Co character (panel d). The reduction of Cu d ↓ occupation and the enhancement of Co d ↓ occupation already hints at a transfer of d ↓ electrons from Cu to Co. This assumption shall be verified further in Section 3.3.
The quite complicated occupation-number profiles reproduced in Figure 1 translate into smooth magnetisation spectra (Figure 3). The ferromagnetic part of the sample ('Co', cyan line) becomes considerable demagnetized to about 50 per cent of its initial value already during and directly after the laser pulse. A fraction of this demagnetisation is caused by spin-conserving transfer of magnetic moment from Co into  Cu (cf. Section 3.3). Hence the demagnetisation of the entire sample (black spectrum, 'total') is weaker and follows the demagnetisation of Co with a time delay.
The magnetisation of the Cu region (magenta in Figure 3) shows a maximum at about 10 fs and likewise indicates transfer of magnetic moment from Co into Cu across the interface. A closer analysis reveals that shortly after the pulse Cu atoms gain magnetic moment which subsequently is lost due to spin-flip processes mediated by the bath or by spin-orbit coupling. We note in passing that the magnetisation dynamics in Co and Cu are very weakly modulated by the laser pulse in the range from about t = −5 fs to about 15 fs.

Spin transfer across the interface
An established picture of spin transport is that majority electrons in the magnetic region of the sample are excited by the laser pulse and then propagate across the interface into the nonmagnetic region [5,6,9,24]. Another mechanism is optically induced inter-site spin transfer (OISTR; Ref. [12]). Here, a laser pulse causes a rearrangement of spin polarisation at the interface which can reverse the magnetic moments of layers. In agreement with Chen et al [25], our simulations establish a similar mechanism: a reflow from Cu d orbitals into empty Co minority-spin d orbitals across the interface. To study this effect we simulated a cluster of 10 Co and 30 Cu atoms with periodicity parallel to the interface. This setup resembles a thin fcc Co layer on a Cu substrate with a (001) interface and facilitates the analysis of currents.
The reflow is illustrated by means of occupation numbers of the interfacial sites. The spin-resolved occupation numbers of Co sp and of Cu sp orbitals are strongly modulated by the laser pulse (Figure 4). A closer look reveals that simultaneously p orbitals become populated, while d orbitals become depopulated, which is readily explained by dipole transitions (δl = −1; both p and d occupation numbers exhibit similar slopes but with opposite sign). The s orbitals become populated slightly after the p orbitals because first p orbitals have to be populated (before the laser pulse both s and p orbitals are weakly populated). The previous effect is contrasted by a depopulation of the d orbitals, except for most of the Co d ↓ orbitals; see panel (b). The occupation numbers of one of the Co and Cu d ↓ orbitals exhibit similar depletion shapes, which is interpreted by a dipole transition into respective p orbitals. The remaining d orbitals behave similarly but with Co populated and Cu depopulated: this observation suggest a spin-conserving 'reflow' from Cu orbitals across the interface into these Co orbitals. Spin-resolved currents, discussed below, support this interpretation.
For the reflow effect it is crucial that the involved orbitals provide a sizable density of states: respective computations for aluminum (instead of copper) exhibit a strongly reduced effect, which is in line with the small density of states of aluminum in the relevant energy range.
The previous considerations corroborate the importance of interfaces for ultrafast spin transport. The local imbalance of occupation at an interface enhances the effect of the laser pulse; an interface may therefore be regarded as an additional 'source of disequilibrium' §. The findings reported so far suggest also to distinguish two induced currents within the sample: a flow of hot, highly mobile s and p electrons contrasted by a reflow of 'tepid', weakly mobile d electrons. These two contributions can be distinguished § We note in passing that simulations for homogeneous systems exhibit considerably less demagnetisation.  The charge currents, defined in Equation (6), are evenly distributed within the sample [ Figure 5(a)]. They are significantly modulated during the laser pulse; confer the alternating blue and red stripes at about 0 fs.
Regardless of the rather uniform distribution of the charge currents across the sample, spin is transferred from the Co region into the Cu region [panel (b)]. The origin of the spin currents is clearly associated with the interface (dark red area around 0 fs), from where it expands covering almost the entire sample. This expansion happens at a high velocity, as indicated by the dashed red line.
The redistribution of d minority-spin electrons at the interface at t ≈ 0 fs causes two phenomena. First, the currents in the Co regions are strongly modulated, which is explained by scattering of electrons at the sample's Co edge and at the interface [cf. the Co region in panel (d) in the range from 5 fs to 40 fs with its intricate pattern]. Second, there is the reflow effect which shows up clearly as elongated blue area in the Cu region. It can be viewed as a depletion region of minority-spin d electrons which, starting at the interface, 'invades' the Cu region with low velocity (dashed blue line). It vanishes after about 40 fs due to dephasing. Moreover, this reflow -figuratively speaking: a domino reaction -contributes weakly to the total spin current shown in panel (b); notice within this respect the different color scales for (b) and (d).

Outlook
In this paper we have delineated that a density-operator approach based on a tightbinding description of the electronic structure is able to capture and reproduce essential phenomena of ultrafast spin dynamics and transport in inhomogeneous systems. To come closer to experiments, however, some areas need improvement. To name but a few: the coupling to the bath could include electron-phonon coupling strengths that are calculated for the actual sample (material) and the tight-binding parameters could be optimized for inhomogeneities.
Concerning extensions of the approach, it is conceivable to incorporate electronelectron scattering and the electrostatic interaction, in particular screening, between regions with in-or decreased population (charge surplus or shortage). Noncollinear spin textures open up additional transport channels (that is spin-flipping hopping in the HamiltonianĤ 0 ) which could cause stronger demagnetisation [25] and stronger spin currents; the textures may be either intrinsic (as in, e. g., Mn 3 Sn [26,27] and Mn 3 Ir [28] or thermal (that is fluctuating magnetic moments, taken for example from atomistic spin dynamics simulations [29]).
In order to deduce spin-resolved currents, we extend the site indices to site, orbital, and spin: k → (k, α, σ) with σ =↑, ↓ gives j kl = i 2 in which Σ µ is a block Pauli matrix and [·, ·] + is the anti-or +-commutator.