Time-Dependent Multiconfiguration Self-Consistent-Field Study on Resonantly Enhanced High-Harmonic Generation from Transition Metal Elements

We theoretically study high-harmonic generation (HHG) from transition metal elements Mn and Mn$^+$, using full-dimensional, all-electron, first-principles simulations. The HHG spectra calculated with the time-dependent complete-active-space self-consistent-field (TD-CASSCF) and occupation-restricted multiple-active-space (TD-ORMAS) methods exhibit a prominent peak at $\sim 50$ eV, successfully reproducing resonant enhancement observed in previous experiments [R.A. Ganeev \textit{et al.}, Opt. Express \textbf{20}, 25239 (2012)]. Artificially freezing $3p$ orbitals in simulations results in its disappearance, which shows the essential role played by $3p$ electrons in the resonant harmonics (RH). Further transition-resolved analysis unambiguously identifies constructively interfering $3p$-$3d$ ($m=0,\pm 1$) giant resonance transitions as the origin of the RH, as also implied by its position in the spectra. Time-frequency analysis indicates that the recolliding electron combines with the parent ion to form the upper state of the transitions. In addition, this study shows that the TD-CASSCF and TD-ORMAS methods can be applied to open-shell atoms with many unpaired inner electrons.


I. INTRODUCTION
High intensity and ultrashort laser pulses have become an indispensable tool both in scientific researches and industrial applications for studying or manipulating the properties of matter. High-harmonic generation (HHG) is one of the important research domains that have emerged thanks to the remarkable advancement of highintensity ultrashort laser technologies. HHG is a nonperturbatively nonlinear optical process in which a fundamental strong-laser field is converted into harmonics of very high orders upon interaction with atoms, molecules, and solids. The nature of the HHG process is closely intertwined with the electronic structure and dynamics of the generating medium. Hence, a variety of quantum scale phenomena have been succesfully identified by devising specialized measurement techniques such as electronic structure detection [1], observation of Rabi flopping [2], multi-channel interference [3], and spectroscopy of Cooper minimum [4,5].
High-harmonic (HH) radiation is an excellent source of coherent XUV photons that can fit into a labroom [6][7][8]. A number of applications for studying atomic and material properties have been reported [1][2][3][4][5] demonstrating its potential as a reliable tool for studying light-matter interaction in the attosecond time scale. Further increase in harmonic intensity is desirable to fully explore its areas of use.
It has frequently been reported that the use of transition metal plasma as a generating medium leads to resonant enhancement of a single or a few harmonics [9][10][11][12][13][14][15][16]. The resonant harmonics (RH) lie close to the giant transition lines in the absorption or emission response of the elements [9], e.g., 27 eV in Sn plasma [10,13,14,17], 20 eV in In [9], and 50 eV in Mn [12]. The RH generation would offer an attractive way to increase HHG yield around a certain photon energy. It might also serve as a new platform to explore multielectron dynamics in intense laser fields using HHG, which is usually considered to be of single electron nature in most cases.
In this work, we theoretically investigate the resonant HHG process, using Mn and its cation Mn + as target systems for the scrutiny of the underlying mechanism. The resonance in HHG spectra from Mn + was observed experimentally by Ganeev et al. [12] where the harmonic peak around 50 eV is enhanced by more than an order of magnitude relative to neighboring harmonics, and the 3p-3d resonance was suggested to be relevant, as implied by the peak position [12].
Prior to the present work, there have been theoretical efforts on RH from transition metal elements over the past decade. Milošević in Ref. [18] and [19] studied the effect of coherent superposition in the initial state and found that a three-step process starting from an excited state but returning to the ground one exhibits an enhanced harmonic. In Ref. [20], the line shape of resonant harmonic is discussed in terms of Fano lineshape. A modelling of the autoionizing state is performed in Refs. [16,17,21,22] using a parametrized potential barrier. Milošević has also investigated the property of resonant harmonic such as intensity and phase in Ref. [23]. Other reports of varying elaboration have also been published, see e. g. Ref. [24][25][26]. Most of the above-mentioned attempts use an effective model potential within the singleactive-electron approximation. Although not targeted at transition metal plasma, one particular work that starts to consider the possibility of multielectron effects has been done by Redkin and Ganeev [27], who have simulated a fullerene-like model system using multiconfiguration time-dependent Hartree-Fock (MCTDHF) method within the two-active-electron jelliumlike sphere approximation.
In the present work, we do all-electron threedimensional (3D) ab initio simulations based on the timedependent multiconfiguration self-consistent-field (TD-MCSCF) methods [28], which describe the system wavefunction by the superposition of Slater determinants consisting of time-dependent spin-orbital functions. Specifically, we apply state-of-the-art implementation of the time-dependent complete-active-space self-consistent-field (TD-CASSCF) [29][30][31] and timedependent occupation-restricted multiple-active-space (TD-ORMAS) [32] methods, which classify spatial orbitals into doubly occupied core and correlated active orbitals. Previously, these methods have been applied to either closed-shell systems or systems having a single unpaired valence electron [30,[32][33][34]. Here, we extend our methods to general open-shell atoms such as transition metals, having many unpaired electrons (5 and 6 for Mn and Mn + , respectively) that can equally participate in the dynamics under strong laser fields. We successfully reproduce the RH at ∼ 50 eV and unambiguously identify the 3p-3d giant resonance as its origin, by taking full advantage of TD-CASSCF and TD-ORMAS to analyze transition dynamics between different orbitals. Our results show that the three 3p-3d lines (m = 0, ±1) constructively interfere to form the RH peak.
This paper is organized in the following way. The overview of the two methods used for our simulations is given in Sec. II. The results are presented and discussed in Sec. III. Conclusions and future possibilities are given in Sec. IV. Hartree atomic units are used throughout unless otherwise noted.

II. TD-CASSCF AND TD-ORMAS METHODS
We consider an N -electron atom (or ion) with atomic number Z irradiated by a laser field E(t) linearly polarized along the z axis. In the velocity gauge and within the dipole approximation, its dynamics is described by the time-dependent Schrödinger equation, a Experimental ionization potential in eV [36]. b Barrier-suppression intensity in W/cm 2 . c Ground state configuration [36].
with A(t) = − E(t)dt being the vector potential.
In the TD-CASSCF and TD-ORMAS methods, we express the total wave function as, whereÂ denotes the antisymmetrization operator, Φ fc and Φ dc the closed-shell determinants formed with n fc time-independent doubly occupied frozen-core and timedependent doubly occupied n dc dynamical-core orbitals, respectively, and {Φ I } the determinants constructed from n a active orbitals. Whereas in the TD-CASSCF method the active electrons are fully correlated among the active orbitals within prescribed numbers of up-and down-spin electrons, the TD-ORMAS method further subdivides the active orbitals into an arbitrary number of subgroups, specifying the minimum and maximum number of electrons accommodated in each subgroup. We specifically consider Mn and Mn + in the present study, whose ionization potential, barrier-suppression intensity [16,35], and the ground-state configuration are summarized in Table I. Orbital subspace decomposition used in this study is shown in Fig. 1. Note that at least 15 spatial orbitals are required for the correct spin multiplicities. TD-CASSCF simulations use 52 and 44 determinants for Mn and Mn + , respectively. In TD-ORMAS simulations, up to two-electron excitations from Active1 to Active2 are allowed, which results in 86510 determinants for Mn and 66068 for Mn + .
The equations of motion (EOMs) describing the temporal evolution of the CI coefficients {C I (t)} and the orbitals {ψ p (t)} are derived on the basis of the timedependent variational principle (TDVP) [37][38][39][40] and read, whereĤ denotes the total Hamiltonian,ĥ the one-body Hamiltonian,Q = 1− q |ψ q ψ q | the projector onto the orthogonal complement of the occupied orbital space.F is a non-local operator describing the contribution from the interelectronic Coulomb interaction, defined aŝ where D and P are the one-and two-electron reduced density matrices, andŴ r s is given, in the coordinate space, by The matrix element R q p is given by, with h q p = ψ q |ĥ|ψ p . R q p 's within one orbital subspace (frozen core, dynamical core and each subdivided active space) can be arbitrary Hermitian matrix elements, and in this paper, they are set to zero. On the other hand, the elements between different orbital subspaces are determined by the TDVP. Their concrete expressions are given in Ref. [32], where iX q p = R q p + h q p is used for working variables.
Our numerical implementation [30] employs a spherical harmonics expansion of orbitals with the radial coordinate discretized by a finite-element discrete variable representation [41][42][43][44]. Specifically, to obtain the ground states we use 12 finite elements each of which contains 25 grid points. For subsequent time-dependent simulations, up to 62 finite elements (including those used in the ground state) are employed. The initial ground state is obtained through imaginary time propagation of the EOMs. The Hartree-Fock energies (−1149.866252 a.u. for Mn and −1149.649383 a.u. for Mn + ) perfectly match the values reported in Ref. [45]. For the CASSCF and ORMAS cases, electron correlation leads to the lowering of the ground-state energy. For example, the OR-MAS method yields the ground state energy of Mn to be −1150.0760 a.u.. We confirm that the total orbital and spin angular momenta of the ground state match the term notations for both Mn ( 6 S 5/2 ) and Mn + ( 7 S 3 ).
We calculate HHG spectra as the magnitude squared of the Fourier transform of the dipole acceleration a(t), defined as [30], Here, the additional term ∆(ṗ z ) accounts for the correction to the Ehrenfest formula in the presence of frozen core orbitals. Its explicit expression is found in Ref. [30].

III. RESULTS AND DISCUSSIONS
A. Resonant high-harmonic emission from Mn and Mn + Let us first examine if our simulations reproduce the resonant harmonics in the HH response of Mn and Mn + . The harmonic spectra from Mn and Mn + obtained with the TD-CASSCF method for a fundamental laser field with 770 nm central wavelength, 3 × 10 14 W/cm 2 peak intensity, and foot-to-foot four-cycle sin 2 pulse shape are shown in Fig. 2(a) and (b) (blue solid, marked as "fz. 2p"), respectively. The results of the MCTDHF simulations, in which all the 15 orbitals in Fig. 1(a) and (b) are treated as active, are also shown (green thick solid curves). The perfect overlap of the results by the two methods indicates numerical convergence for this number of orbitals. The results of the TD-ORMAS simulations are plotted in Fig. 3 for the same laser parameters as in Fig. 2. In both Figs. 2 and 3, we clearly see an RH slightly above 50 eV, substantially enhanced in comparison with neighboring harmonics, for both Mn and Mn + , whose position is in excellent agreement with the experimental value (∼ 50 eV [12]).
We have calculated HH spectra for 1333 nm wavelength and 10 14 W/cm 2 peak intensity with the TD-ORMAS method (Fig. 4), while keeping the ponderomotive energy unchanged. In spite of substantial difference in laser parameters, we can see that the resonant harmonic peak remains at ∼ 50 eV. Thus, the RH position is governed by atomic properties, rather than laser parameters. Then, we make use of the flexibility in orbital-subspace decomposition to find out which orbitals contribute to RH. We have performed TD-CASSCF simulations by varying the boundary between the active and frozen spaces in Fig. 1. Whereas freezing 3s virtually does not change the spectrum, freezing up to 3p leads to the disappearance of the RH peak (Fig. 2). This indicates that the appearance of the enhanced peak involves the dynamics of 3p electrons.
In Fig. 5 are shown the (single-photon) excitation spectra of Mn, Mn + , and Mn 2+ , obtained as a Fourier transform of the dipole response to a quasi-delta-function pulse with the field being finite at three time steps (10 8 W/cm 2 peak intensity). Although there are slight differences between the TD-CASSCF and TD-ORMAS results, we see a strong excitation line at ∼ 50 eV in all the cases, which reproduces the position of the well-known 3p-3d giant resonance line [48,49] and coincides with the RH in the HHG spectra. These observations strongly suggest that the RH originates from the 3p-3d resonance line, as also implied experimentally [12].
Before ending this subsection, let us briefly discuss the cutoff energies. The arrows in Fig. 3 mark the cutoff positions E c expected from the cutoff law E c = I p + 3.17U p with U p being the ponderomotive energy. Even if the simulation starts from Mn or Mn + , the HHG spectra extend further up to the cutoff corresponding to Mn 2+ . In- deed, as expected from the barrier-suppression intensity (Table I) and confirmed by Fig. 6 showing the temporal variation of the fraction of each species, Mn is mostly ionized in the early stage, and Mn + is further substantially ionized to Mn 2+ . The comparison between Fig. 6 and the time-freqency structure of HHG shown in Fig. 7 also indicates that the higher plateau appears after the production of Mn 2+ . In spite of its high ionization potential, harmonic response of Mn 2+ is enhanced probably through laser-induced electron recollision [33,34,50].

B. Transition-resolved analysis
The results in the previous Subsection motivate us to analyze contributions from individual transition lines. For the transition-resolved analysis, let us rewrite the dipole acceleration Eq. (8) as, wheref = −Z(z/r 3 ), and {|n } denotes the initial orbitals, obtained through imaginary-time relaxation. Since each initial orbital has a definite parity, the terms for m = n vanish. Then, we can view, as a contribution from a transition between orbital pair m and n to the dipole acceleration.
With the orbitals used in TD-ORMAS simulations [ Fig. 1(c) and (d)], we can identify the following 18 transitions satisfying the selection rule: where the subscripts denote the magnetic quantum numbers. We have calculated the power spectrum of each transition line as the magnitude squared of Fourier transform of α(m, n, t).
The spectrum of each of 3p 0,±1 ↔ 3d 0,±1 , the contribution of their sum, and the total spectrum of all the lines listed in Eq.(11) are presented in Fig. 8 for Mn and Mn + . The sum contribution of 3p 0,±1 ↔ 3d 0,±1 has a clear peak at ∼ 50 eV and dominates the total contribution in that photon energy region. The other transitions involving 3p such as 3p 0 ↔ 4s are by orders of magnitude weaker. Moreover, the position and form of the peak agree very well with those of the RH in the HHG spectra shown in Fig. 3. It should also be noticed that the sum spectrum of the three 3p 0,±1 ↔ 3d 0,±1 lines is approximately one order of magnitude stronger than the contribution of each transition. These observations unambiguously establish that the resonant harmonic at ∼ 50 eV, experimentally discovered [12] and numerically reproduced by our ab initio simulations is driven by a constructive interference of electron dynamics occurring between 3p and 3d orbitals.
Superposed on the spectrograms in Fig. 7, we plot, as a function of time, the recombination energy, defined as the sum of the kinetic energy of the returning electron and the ionization potential, or, harmonic photon energy expected from the three-step model [46,47]. We see that RH photons are emitted mainly when the recombination energy is ∼ 50 eV, where the recolliding electron has to be recaptured by the parent ion in order to induce a 3p → 3d transition. This is consistent with recombination to autoionizing states (the upper states of the 3p-3d giant resonance lines in the present case), proposed in Refs. [16,17,21,22]. Whereas these studies used the single-activeelectron approximation with a model potential barrier, our all-electron ab initio simulations support this process as a mechanism of RH generation from Mn and Mn + . 3p -1 <--> 3d -1 3p 0 <--> 3d 0 3p +1 <--> 3d +1 sum of 3p <--> 3d sum of all 3p 0 <--> 4s 3p -1 <--> 3d -1 3p 0 <--> 3d 0 3p +1 <--> 3d +1 sum of 3p <--> 3d sum of all 3p 0 <--> 4s on the experimentally observed resonance harmonics [12]. While its position suggests the relevance with the 3p-3d giant resonance lines, we have performed a series of analyses to unambiguously verify it. First, we have taken advantage of flexibility in orbital subspace decomposition to vary the boundary between frozen-core and active orbitals, and found that freezing up to 3p leads to the disappearance of the RH, which shows the essential role of the 3p electrons. Then, we have calculated the contribution of each transition between initial orbitals to harmonic spectra. It has indeed revealed that the RH is dominated by the 3p-3d (m = 0, ±1) transitions, constructively interfering. It has followed from the inspection of the HHG time-frequency structure that the RH is emitted mainly, if not exclusively, when the sum of the kinetic energy of the returning electron and the ionization potential of Mn or Mn + is ∼ 50 eV. This implies that the electron recombines to the autoionizing upper state of the 3p-3d transitions, as proposed previously within the single-active-electron approximation using a model potential barrier [16,17,21,22].