text:
# Optical creation of vibrational intrinsic localized modes in anharmonic lattices with realistic interatomic potentials

## I Introduction

Over the past several years, theoretical studies of the dynamics of anharmonic periodic lattices have established the existence of intriguing vibrational excitations, characterized by well-localized displacement patterns. These so-called intrinsic localized modes (ILMs) can exist at any site in a perfect lattice, in contrast to localized impurity modes in harmonic defect crystals. While vibrational ILMs with atomic scale localization have been obtained for increasingly realistic lattice-dynamical models, their experimental study has been impeded by the lack of direct methods for their excitation and verification. Recently, ILMs in a complex quasi-1D charge-density wave system were inferred from resonance Raman data, through the use of a coupled electron-vibration model restricted to a single repeat unit.

Previously we have demonstrated theoretically that driven ILMs can exist as a steady-state response to an applied spatially homogeneous monochromatic driving force. As a natural extension of that work, and to address a key experimental question, we here describe theoretically an avenue for the transient optical creation of ILMs: we show how they can be produced in a 1D model lattice with realistic potentials by means of laser pulses whose time dependence is designed by an efficient optimal control scheme.

Owing to their high power densities, lasers are attractive sources for exciting large-amplitude ILMs. Indeed, there are some notable examples of experiments with powerful lasers in the regime of anharmonic lattice vibrations. For instance, phonon resonances measured in the ferroelectric LiNbO3 by experiments using single visible laser pulses with a duration of 60 fs and an energy per pulse of 5 $`\mu `$J were interpreted in terms of “overtones” of the very anharmonic lowest-energy TO phonon of $`A_1`$ symmetry. Also, experiments on Ti2O3 using single visible pulses with 10 nJ energy and 70 fs widths to excite the $`A_{1g}`$ mode apparently resulted in atomic displacements of about 0.07 Å – corresponding to 2% of the interatomic spacing – and revealed anharmonic behavior.

Of more direct relevance here is the fact that experimental laser pulse shaping techniques  provide considerable flexibility in the time dependence of the applied force. Indeed, the use of tailored fields for vibrational excitation has attracted much interest, mainly in the context of optical control of dissociation and reactions in molecular chemistry, but also for the selective excitation of optical phonons in time-domain spectroscopy. We will focus on two methods for the transient optical creation of ILMs: (i) impulsive stimulated Raman scattering (ISRS) excitation  by a sequence of femtosecond pulses at THz repetition rates from a laser operating at visible or near-visible frequencies, and (ii) infrared (IR) excitation by a picosecond far-IR laser pulse. For both mechanisms, the system’s dynamical response can be described classically, provided the underlying laser frequency for the ISRS case is well off resonance with vibrational and electronic transitions.

Since ILMs are complex, large-amplitude vibrational excitations, the external fields necessary for their creation are likely to have a complicated time dependence and large magnitudes. It is therefore advantageous to determine the optimal external fields by a systematic scheme. In engineering, the analogous task of designing the time dependence of an external force to steer a dynamical system towards a desired target state is a fundamental problem. Optimal control theory  provides a solution with a rigorous mathematical foundation, based on the variational minimization of a positive objective functional. In the realm of atomic dynamics, this approach has been successful in the design of external electric fields for selective bond excitation in models of small harmonic lattices  and small anharmonic molecules. We apply a similar scheme to the creation of ILMs.

The following section discusses our theoretical framework by providing the necessary details of the optimal control scheme and describing the specific anharmonic model lattice we use. Our numerical results concerning the optical creation of ILMs are given in Sec. III, which is divided into two parts according to the two different optical methods we consider. In Sec. IV we address aspects of the experimental feasibility and discuss our results. Section V concludes the paper, and an appendix provides additional qualitative insight. Some of the results presented here were summarized in a letter.

## II Theoretical background

Although some notable exceptions have appeared, most studies of ILMs have considered 1D or 2D model lattices. This restriction to simpler model systems has facilitated progress towards a theoretical understanding of basic ILM properties, without the additional numerical complications encountered with more realistic 3D models of crystals. In addition to the reasonable expectation that many of these basic properties will transfer to the 3D case, there is evidence that 1D models may apply directly to some types of motion in 3D crystals. For instance, in Ref.  a realistic 3D lattice-dynamical model was considered and ILMs were obtained with displacement patterns localized along the edge of the crystal, which can be regarded as a natural generalization of a 1D lattice. Furthermore, it is well known that along some high symmetry directions in 3D crystals, the harmonic lattice dynamics map onto an effective 1D model involving the collective motion of lattice planes. We have performed preliminary studies showing that this mapping also occurs in the anharmonic case, for certain polarization directions. This point is addressed in more detail later. Hence, for a demonstration of the optical creation of ILMs, we consider a diatomic 1D model lattice that incorporates realistic features, such as standard interparticle potentials and measured harmonic properties of real crystals.

### A The system: Hamiltonian and dynamics

For longitudinal motion in an externally driven 1D system, the Hamiltonian is

$$H=\underset{n}{}\left[\frac{p_n^2}{2m_n}+\underset{l>0}{}V_{n,nl}\left(r_nr_{nl}\right)f_n^{\mathrm{ext}}(t)r_n\right],$$
(1)
where particle $`n`$ has mass $`m_n`$, position $`r_n`$ and momentum $`p_n`$ and interacts with particle $`nl`$ via a potential $`V_{n,nl}(r)`$, to be specified below. The external force is given by $`f_n^{\mathrm{ext}}(t)=\frac{1}{2}𝒫_n^2(t)`$ and $`f_n^{\mathrm{ext}}(t)=q_n(t)`$ for ISRS and IR excitation, respectively. Here $`(t)`$ is the longitudinally polarized electric field, $`𝒫_n(𝒫/r_n)_0`$ is the electronic polarizability derivative evaluated at the equilibrium configuration, and $`q_n`$ is the effective charge. With $`𝒫_n=(1)^n𝒫`$ and $`q_n=(1)^nq`$, the external forces for both excitation methods have equal magnitudes and alternating signs:

$$f_n^{\mathrm{ext}}(t)=(1)^n(t).$$
(2)
The vibrational dynamics of this 1D system are described by Hamilton’s equations

$$\dot{r}_n=\frac{H}{p_n}=\frac{p_n}{m_n}$$
(4)
$$\dot{p}_n=\frac{H}{r_n}=f_n(\{r_m\})+(1)^n(t),$$
(5)
where

$`f_n(\{r_m\}){\displaystyle \underset{l>0}{}}`$ $`[V_{n,nl}^{}(r_nr_{nl})`$ (7)
$`V_{n+l,n}^{}(r_{n+l}r_n)]`$
is the total internal force on particle $`n`$, with $`V_{n,m}^{}(r_nr_m)(dV_{n,m}/dr)|_{r=r_nr_m}`$. From Eqs. (2) we readily obtain the equations of motion

$$m_n\ddot{r}_n=f_n(\{r_m\})+(1)^n(t),$$
(8)
which are more convenient for numerical simulations of the system dynamics.

### B The scheme: optimal control theory

In the following, we describe the most important aspects of the optimal control scheme used in this work. More details can be found in Refs.  and , whose compact notation we adopt with some modifications. We consider a system of $`N`$ particles. Unless explicitly stated otherwise, $`N`$-dimensional vectors are denoted by lower case bold Roman letters and $`N\times N`$-dimensional matrices by upper case bold Roman letters. Moreover in phase space, $`2N`$-dimensional vectors are denoted by lower case bold Greek letters and $`2N`$$`\times `$$`2N`$-dimensional matrices by upper case bold Greek letters.

We combine the $`N`$-dimensional vectors $`𝐫`$ and $`𝐩`$ for the particle positions and momenta, respectively, to form a $`2N`$-dimensional phase space vector

$$𝝃^T(𝐫^T,𝐩^T)(r_1,\mathrm{},r_N,p_1,\mathrm{},p_N),$$
(9)
where the superscript $`T`$ denotes the transpose. Equations (2) are then rewritten as

$$\dot{𝝃}=\mathit{\varphi }[𝝃,(t)],$$
(10)
with

$$\mathit{\varphi }^T[𝝃,(t)][(𝐌^1𝐩)^T,𝐟^T+\stackrel{~}{𝐪}^T(t)],$$
(11)
where $`𝐌`$ is the $`N`$$`\times `$$`N`$-dimensional diagonal mass matrix with elements $`(𝐌)_{nl}m_n\delta _{nl}`$, $`(𝐟)_nf_n(\{r_m\})`$, and the definition of the coupling vector $`\stackrel{~}{𝐪}`$ with components $`(\stackrel{~}{𝐪})_n(1)^n`$ allows a compact description of the external force terms. With this notation, we can rewrite Eq. (8) as

$$𝐌\ddot{𝐫}=𝐟+\stackrel{~}{𝐪}(t).$$
(12)
We furthermore specify initial conditions

$$𝐫(t=0)=𝐫_i,$$
(14)
$$𝐩(t=0)=𝐩_i$$
(15)
at $`t=0`$.

In order to apply optimal control methods to our problem, we first need to define a positive functional that reflects the physical objectives to be reached. Starting from the lattice in some initial configuration $`(𝐫_i,𝐩_i)`$, we want to excite a given anharmonic mode at a specified final time $`t_f`$, while keeping the magnitude of the external force within reasonable limits. Combining the target mode positions $`𝐫_f`$ and momenta $`𝐩_f`$ into a final phase space vector $`𝝃_f`$, we define the objective functional

$`J[𝝃,(t)]=`$ $`{\displaystyle \frac{1}{2}}\left[𝝃(t_f)𝝃_f\right]^T𝚿\left[𝝃(t_f)𝝃_f\right]+`$ (17)
$`{\displaystyle \frac{1}{2}}\psi _{}{\displaystyle _0^{t_f}}𝑑t^2(t),`$
where the nonzero elements $`(𝚿)_{\alpha \alpha }\psi _\alpha `$ of the $`2N`$$`\times `$$`2N`$-dimensional diagonal matrix $`𝚿`$ are positive weight factors, as is $`\psi _{}`$. Note that in order for all terms in the objective functional to have the same units, not all of the weight factors can be unitless. Furthermore, only trajectories $`𝝃(t)`$ satisfying the equations of motion (10) are admissible during the optimization. Including this as a constraint in the objective functional, we obtain the modified objective functional

$$\overline{J}[𝝃,(t)]=J[𝝃,(t)]_0^{t_f}𝑑t𝝀^T\left\{\dot{𝝃}\mathit{\varphi }[𝝃,(t)]\right\},$$
(18)
where $`𝝀^T=(\lambda _1,\mathrm{},\lambda _{2N})`$ is a $`2N`$-dimensional vector of time-dependent Lagrange multipliers. This modified objective functional is minimized with respect to the trajectories $`𝝃`$ and external force $`F(t)`$ to obtain the optimal force. For clarity, these quantities are treated separately in the following two paragraphs.

Variational minimization with respect to the trajectories $`\{\xi _\alpha (t)\}`$, including an integration by parts of the term $`_0^{t_f}𝑑t𝝀^T\delta \dot{𝝃}`$, yields dynamical equations for the Lagrange multipliers

$$\dot{𝝀}+𝚽𝝀=0,$$
(19)
where $`𝚽`$ is a $`2N`$$`\times `$$`2N`$-dimensional time-dependent matrix with elements $`(𝚽)_{\alpha \beta }\varphi _\beta [𝝃,(t)]/\xi _\alpha `$. At the final time $`t_f`$, these dynamical equations are subject to boundary conditions

$$𝝀(t_f)=𝚿[𝝃(t_f)𝝃_f].$$
(20)
For our definition of $`\mathit{\varphi }[𝝃,(t)]`$ \[see Eq. (11)\], the matrix $`𝚽`$ decomposes into four $`N`$$`\times `$$`N`$-dimensional blocks:

$$𝚽=\left(\begin{array}{cc}\mathrm{𝟎}& 𝐊\\ & \\ 𝐌^1& \mathrm{𝟎}\end{array}\right),$$
(21)
where we have defined the $`N`$$`\times `$$`N`$-dimensional time-dependent, symmetric dynamical matrix $`𝐊`$ with elements $`\left(𝐊\right)_{nl}f_n(\{r_m\})/r_l`$. We can now simplify the description by considering separately the Lagrange multipliers $`(𝝀^r)^T(\lambda _1,\mathrm{},\lambda _N)`$ for the positions and $`(𝝀^p)^T(\lambda _{N+1},\mathrm{},\lambda _{2N})`$ for the momenta. Then, Eq. (19) becomes

$$\dot{𝝀}^r𝐊𝝀^p=0,$$
(23)
$$\dot{𝝀}^p+𝐌^1𝝀^r=0,$$
(24)
which can be combined to yield

$$𝐌\ddot{𝝀}^p=𝐊𝝀^p.$$
(25)
We note that these dynamical equations for the momentum Lagrange multipliers correspond to the equations of motion for a lattice of $`N`$ fictitious particles with masses $`m_n`$ and instantaneous positions $`\lambda _n^p`$, interacting via time-dependent harmonic forces. Combining Eqs. (20) and (21), we obtain the corresponding boundary conditions

$$𝝀^p(t_f)=𝚿^p[𝐩(t_f)𝐩_f],$$
(27)
$$\dot{𝝀}^p(t_f)=𝐌^1𝚿^r[𝐫(t_f)𝐫_f],$$
(28)
where $`𝚿^p`$ and $`𝚿^r`$ are $`N`$$`\times `$$`N`$-dimensional diagonal weight factor matrices with elements $`(𝚿^p)_{nl}=\delta _{nl}\psi _{n+N}`$ and $`(𝚿^r)_{nl}=\delta _{nl}\psi _n`$. Note, that these boundary conditions for the Lagrange multipliers are specified at the final time $`t_f`$. Thus their determination requires the knowledge of the final positions $`𝐫(t_f)`$ and momenta $`𝐩(t_f)`$ of the particles in the actual lattice. We emphasize that the introduction of Lagrange multipliers only serves the purpose of including the dynamics of the driven lattice as a constraint in the objective functional (17). This results in coupled equations for the dynamics of the actual lattice and that of the Lagrange multipliers.

The optimal external force is now obtained by minimizing the modified objective functional (18) with respect to $`(t)`$. This minimization can be done for an external force $`(t)`$ whose time dependence is allowed to be arbitrary  or which has a prescribed analytic form. In view of the important aspect of experimental feasibility, we discuss the latter approach. For ISRS excitation, the laser frequency is neglected and we constrain $`(t)`$ to be a sequence of Gaussian pulses

$$^{\mathrm{ISRS}}(t)=\underset{i}{}S_ie^{(tt_i)^2/\mathrm{\Delta }^2}$$
(29)
with individual heights {$`S_i`$} and pulse center times {$`t_i`$}, and common width $`\mathrm{\Delta }`$. For IR excitation, $`(t)`$ is taken to have a linearly chirped IR frequency under a single Gaussian envelope with fixed width:

$$^{\mathrm{IR}}(t)=Se^{(tt_0)^2/\mathrm{\Delta }^2}\mathrm{sin}(\theta +\omega t+\alpha t^2).$$
(30)
In both cases, the external force depends on a set of variable parameters {$`\tau _j`$}, namely $`\{S_i\},\{t_i\}`$ and $`\mathrm{\Delta }`$ for the ISRS case, and $`S,t_0,\theta ,\omega `$, and $`\alpha `$ for the IR case.

The gradient of the modified objective functional with respect to the external force parameters {$`\tau _j`$} has components

$`{\displaystyle \frac{\overline{J}[𝝃,(t,\{\tau _j\})]}{\tau _i}}=`$ (31)
$`{\displaystyle _0^{t_f}}𝑑t\left[\psi _{}(t,\{\tau _j\})+\left(𝝀^p\right)^T\stackrel{~}{𝐪}\right]{\displaystyle \frac{(t,\{\tau _j\})}{\tau _i}},`$ (32)
where Eq. (11) has been used to write $`𝝀^T\mathit{\varphi }/F(t)=(𝝀^p)^T\stackrel{~}{𝐪}`$. The optimal control force $`^{\mathrm{opt}}(t)`$ is obtained by finding the zero of this gradient. Following Ref. , we use an iterative approach for this nontrivial numerical problem. Using an educated guess for the force parameters {$`\tau _j`$}, we integrate the dynamical equations (12) for the actual particles’ driven motion forward in time from $`t=0`$ to $`t=t_f`$, starting from the initial conditions (12). This yields $`𝐊(t)`$ for $`t`$ in the interval $`[0,t_f]`$, as well as the final positions $`𝐫(t_f)`$ and momenta $`𝐩(t_f)`$, which are used to evaluate the boundary conditions (25) for the Lagrange multipliers at $`t=t_f`$. We next integrate Eq. (25) backward in time from $`t=t_f`$ to $`t=0`$. From this we obtain $`𝝀^p(t)`$ for $`t`$ in $`[0,t_f]`$, which is then used to evaluate the gradient components (31). We adapted a fifth-order Gear predictor-corrector molecular dynamics (MD) method  for the time evolution of the particles in the actual lattice and for the Lagrange multipliers. The force parameters are updated within a conjugate gradient scheme, and the procedure is repeated until the gradient (31) is zero to within a specified tolerance.

For the case of purely harmonic potentials $`V_{n,nl}(r_nr_{nl})`$, $`𝚽`$ and $`𝐊`$ are time-independent matrices. Consequently, for initial conditions $`𝐫_i=0`$ and $`𝐩_i=0`$ corresponding to the lattice at rest in its equilibrium configuration, the exact optimal control force can be obtained by direct algebraic manipulation of the matrix equations (10) and (21), without iteration. This is detailed in Refs.  and  for the application to harmonic molecules. Within this harmonic limit, we have used this algebraic approach as an independent test of our iterative MD method described above.

The harmonic limit also allows considerations of the controllability of the system. A system is said to be controllable if an arbitrary specified target state $`𝝃_𝒇`$ can be reached exactly with a suitable external force $`(t)`$. In Ref.  it was pointed out that a harmonic system is completely controllable only if all normal modes couple to the external force. Still, even if a system is not controllable in this rigorous sense, i.e., the external force couples only to a restricted set of normal modes, it may be possible to reach a final state $`𝝃(t_f)`$ that is close, although not exactly equal, to the specified target state $`𝝃_f`$. We will return to this aspect in Sec. III below.

### C The model: interatomic potentials and characteristics

We consider a 1D diatomic lattice with masses $`m`$ and $`M`$ ($`>m`$), where nearest neighbors interact via Born-Mayer plus Coulomb (BMC) potentials

$$V_{mM}(r)=\lambda _{mM}e^{r/\rho }\frac{q^2}{r}$$
(34)
$$V_{Mm}(r)=\lambda _{Mm}e^{r/\rho }\frac{q^2}{r},$$
(35)
while second neighbors interact via pure Coulomb potentials

$$V_{mm}(r)=V_{MM}(r)=\frac{q^2}{r}$$
(36)
and more distant neighbors are assumed to be non-interacting. Note that the interaction between an atom and its nearest neighbors distinguishes between “left” and “right” neighbors. Accordingly, minimization of the total potential energy of the static lattice leads to asymmetric nearest-neighbor equilibrium separations $`R_0^{mM}`$ and $`R_0^{Mm}(R_0^{mM})`$. Although this asymmetry might appear unusual at first sight, it indeed correctly represents the situation encountered by mapping the collective motion of $`(111)`$ planes in a zinc blende structure crystal onto an effective 1D model lattice. Furthermore, this asymmetry is in fact necessary to properly represent a lattice with first-order Raman active vibrational modes.

In a detailed study, we noted that ILMs should be classified according to their “associated” extended lattice mode (ExM), into which they spatially broaden with decreasing amplitude. Since we are focusing on the optical creation of ILMs, we study a model that exhibits ILMs associated with the optically-active ExM. This intuitively reasonable choice will later turn out to be essential. With an external force such as given by Eq. (2), the anharmonic version of the optical zone-center mode (OZCM) with displacement pattern $`A(\mathrm{},1,m/M,1,m/M,\mathrm{})`$ is both first-order Raman and IR active, where for the former the asymmetry $`R_0^{mM}R_0^{Mm}`$ is necessary, as noted above. We therefore need to find model parameters such that ILMs associated with the OZCM exist. In Ref. , we obtained an ILM existence criterion based on the interplay between fundamental harmonic and anharmonic dynamical properties of the ILM’s associated ExM. For our case of interactions via realistic BMC potentials with their dominant soft anharmonicity, this criterion predicts the existence of ILMs associated with the OZCM in diatomic lattices for which the optical branch of the harmonic dispersion relation has a minimum at $`k=0`$. In order to satisfy this constraint on the harmonic properties of our 1D model lattice, while at the same time ensuring reasonably realistic interatomic forces, we determined our BMC potential parameters by fitting the harmonic dispersion to branches of the measured dispersion for a real crystal along a high-symmetry direction. Since we consider both Raman and IR excitation, we focused on real diatomic crystals with the zinc blende structure, for which the $`𝐤=0`$ transverse optical (TO) phonon is both first-order Raman and IR active. In particular, the measured TO phonon branch of ZnS along the $`<`$111$`>`$ direction exhibits the required minimum at $`𝐤=0`$ (see Ref. ). Using ZnS masses $`m=32.1`$amu and $`M=65.4`$amu, together with $`\lambda _{mM}=3.45\times 10^2`$eV, $`\lambda _{Mm}=2.73\times 10^3`$eV, $`\rho =0.279`$Å, and $`q=0.9e`$, we can approximately match those curves, as shown in the upper left panel of Fig. 1. In the resulting model, $`R_0^{mM}=1.67`$Å and $`R_0^{Mm}=2.95`$ Å, and the maximum of the harmonic phonon gap occurs at the harmonic OZCM frequency $`\omega _0=5.22\times 10^2`$rad/fs.

The optimal control scheme described in Sec. II B requires the positions and momenta of a specified target state. In order to obtain accurate predictions for the stationary solutions to Eq. (8) without external driving, we use the well-established rotating wave approximation (RWA) for the particles’ time dependence, but generalized to include static and second-harmonic contributions as well as oscillation at a mode’s fundamental frequency $`\omega `$. The motion of particle $`n`$ is assumed to be of the form

$$r_n(t)=b_n+c_n\mathrm{cos}(\omega t)+d_n\mathrm{cos}(2\omega t)+r_n^0.$$
(37)
After inserting this ansatz into Eq. (8), we multiply the resulting equations by either unity, $`\mathrm{cos}(\omega t)`$, or $`\mathrm{cos}(2\omega t)`$, and average over a single period. For an $`N`$-particle lattice, this yields a system of $`3N`$ coupled nonlinear equations for the static displacements {$`b_n`$}, the fundamental dynamic displacements {$`c_n`$}, and the second harmonic dynamic displacements {$`d_n`$}:

$$0=\frac{1}{2\pi }_0^{2\pi }𝑑\varphi f_n(\{r_m\}),$$
(39)
$$0=m_n\omega ^2c_n+\frac{1}{\pi }_0^{2\pi }𝑑\varphi \mathrm{cos}\varphi f_n(\{r_m\}),$$
(40)
$$0=m_n(2\omega )^2d_n+\frac{1}{\pi }_0^{2\pi }𝑑\varphi \mathrm{cos}(2\varphi )f_n(\{r_m\}),$$
(41)
where $`\varphi =\omega t`$ and $`f_n(\{r_m\})`$ is defined in (7). Once the boundary conditions are specified, these equations can be solved using standard numerical routines; the solutions in conjunction with Eq. (37) constitute the RWA. We verify our RWA predictions by performing direct MD simulations of Eq. (8), using a fifth-order Gear predictor-corrector method. Traditionally, Born-von Karman periodic boundary conditions are employed for the description of bulk properties of “infinite” lattices. These boundary conditions are implemented by setting $`r_{n+N}r_n+L`$, where $`L`$ is the static-lattice equilibrium length of the supercell. Following the nomenclature of Ref. , we denote them standard periodic boundary condition (StdPBCs). While StdPBCs are convenient to describe infinite periodic lattices, we will also consider finite systems with free-end boundary conditions (FBCs).

Applying the RWA to our model lattice, we find that ILMs associated with the OZCM exist, as predicted by the criterion of Ref. . The ILM frequencies are in the harmonic phonon gap. The middle panel of Fig. 1 gives the RWA static and fundamental dynamic displacements of such an OZCM-ILM at $`\omega =0.97\omega _0`$ for a 22-particle system with FBCs. The displacements are normalized to the average equilibrium separation $`R_0(R_0^{mM}+R_0^{Mm})/2`$. For clarity, the corresponding second harmonic dynamic displacements, whose magnitude is less than 8% of the fundamental dynamic displacements, are shown separately in the bottom panel. That this ILM is associated with a first-order Raman active ExM is reflected by the fact that its fundamental dynamic displacement pattern exhibits no inversion symmetry. The static displacements seen at the ends of this finite system are the result of “surface” relaxation. The upper right panel shows the RWA frequency vs amplitude curves for the ILM and for the OZCM in the corresponding 40-particle lattice with StdPBCs. For each of these two modes the amplitude is given by the magnitude of the largest fundamental dynamic displacement $`c_n`$. Also shown are the results of MD measurements of the ILM and OZCM frequencies, which agree to within less than 0.5% with the RWA predictions. The $`\omega (A)`$ curve for this ILM is indistinguishable from the corresponding curve for the 22-particle FBC lattice of the middle panel, as one would expect from the mode’s high localization.

For our purpose of optical ILM excitation, the dynamical stability properties of the optically-active ExM are important, as will be seen later. We have therefore examined the stability of the OZCM in our model lattice within an RWA-based approach detailed in Refs.  and , but generalized to include perturbations of the second harmonic contributions. This stability analysis assumes infinitesimal displacement and velocity perturbations having an exponential time dependence $`\mathrm{exp}(\lambda t)`$, leading to a linear eigenvalue problem for the growth rates {$`\lambda `$}. In previous studies, we have shown that a dynamical ExM instability with a purely real growth rate $`\lambda `$ is intimately connected with the existence of ILMs associated with the ExM. The OZCM in our model lattice exhibits such an ILM-related instability. In the top panel of Fig. 2, we show the predicted maximum growth rate as a function of the OZCM amplitude. MD measurements based on the “projection method” of Ref.  agree to within 8% with the RWA predictions, as indicated by the diamonds. By decomposing the instability perturbation into its spatial Fourier components, we can extract the wave vector $`(k_p)_{\mathrm{max}}`$ of the fastest-growing component, as discussed in Refs.  and . The bottom panel of Fig. 2 plots $`(k_p)_{\mathrm{max}}`$ as a function of the OZCM amplitude. At zero amplitude, $`(k_p)_{\mathrm{max}}`$ vanishes, and it increases to its maximum allowed value $`\pi /(2R_0)`$ over a restricted range of amplitudes. As discussed in Refs.  and , the corresponding wavelength $`2\pi /(k_p)_{\mathrm{max}}`$ introduces a preferred instability length scale at each amplitude. As the OZCM amplitude is increased from zero, the instability length scale decreases from infinity, reaching a maximum value of $`4R_0`$, after which it remains constant for increasing amplitude. Finite-time MD simulations of unstable OZCMs seeded with the fastest growing instability perturbation for various amplitudes reveal that the instability leads to a breakup of the OZCM into a periodic array of localized ILM-like excitations whose spacing is very close to the preferred instability length. The stability properties of the OZCM will be used in Sec. III B below.

## III Optical creation of ILMs

Before we apply the optimal control scheme of Sec. II B to the optical creation of OZCM-ILMs in the model lattice described in Sec. II C, we consider the question of controllability, which was briefly addressed in Sec. II B. We find that this aspect of our problem is different for different choices of the boundary conditions.

For a harmonic diatomic lattice with StdPBCs, it is well known that the harmonic version of the OZCM is the only normal mode that couples to an optical-like force with a spatial dependence such as given in Eq. (2). Hence, independent of the time-dependence of $`(t)`$, this external force can only excite the OZCM. On the other hand, if we start from zero initial conditions in an anharmonic lattice with StdPBCs, this argument still applies at short times, since for the initially small amplitudes the interactions are dominated by the harmonic terms. It remains true even at longer times for our StdPBC model lattice, because the large-amplitude anharmonic OZCM continues to be an exact stationary solution.

On the other hand, when a finite harmonic system with FBCs is considered, the external force (2) couples with appreciable strength to a set of normal modes. However, the size of this set decreases with increasing size of the system, and in the limit $`N\mathrm{}`$ of an infinite system it is again just the FBC normal mode corresponding to the OZCM in a StdPBC lattice which can be excited optically. From these considerations of the harmonic case, we expect the external force (2) to have better control over a system with FBCs than over a lattice with StdPBCs, when anharmonicity is included. However, with FBCs we also expect to see a dependence on the system size, with increased controllability for smaller systems.

In Sec. III A we show how ILMs can be excited “directly” in a finite system with FBCs. As expected, this approach depends on the size of the system. For the creation of ILMs in crystal lattices, we exploit the fact that the OZCM is unstable and breaks up into ILM-like localized vibrations. Hence, although we can only excite the StdPBC OZCM directly, its decay can produce ILM-like vibrations “indirectly.” This is detailed in Sec. III B.

### A “Direct” excitation of ILMs in finite systems

#### 1 ISRS excitation

We first demonstrate the power of the optimal control scheme by showing our results for ISRS excitation of an OZCM-ILM in a 22-particle system with FBCs. Preliminary studies with simpler model systems revealed that the control scheme becomes more successful for longer control periods and larger numbers of pulses in the applied sequence. Keeping the question of experimentally feasible laser pulse sequences in mind, we choose a control interval $`[0,t_f=50T_0]`$, where $`T_02\pi /\omega _0`$ is the period of the harmonic OZCM. At $`t=0`$ the particles are at rest at their equilibrium positions. The system is then driven by a sequence of 49 Gaussian pulses. The MD time step used during the optimization for this case and for all of the simulations in this paper is $`T_0/100`$. For the target state at $`t_f`$ we specify the RWA-predicted positions and momenta of an ILM at frequency $`\omega =0.97\omega _0`$. The displacement pattern is given in the lower two panels of Fig. 1. Our initial studies with simpler systems also showed that a combination of weight factors $`\psi _n=16(m_n\omega _0/2)`$ and $`\psi _{n+N}=16/(m_n\omega _0/2)`$ for the target positions and momenta of particle $`n`$, respectively, and $`\psi _{}=1`$, guarantees a good balance of the various terms in the objective functional. Using these values in our control algorithm yields the pulse sequence given in the top panel of Fig. 3. The common full width at half maximum (FWHM) $`2\sqrt{\mathrm{ln}2}\mathrm{\Delta }`$ of the pulses is 18 fs, and the amplitudes {$`S_i`$} range from zero to 0.13 eV/Å. The bottom panel shows that the application of this rather complex sequence of pulses in MD produces a strikingly “simple” result, namely the creation of a highly-localized excitation which persists almost unchanged well after the applied field ends at $`t=t_f`$.

We note that the pulse sequence shown in the top panel of Fig. 6 does not represent a unique solution of the minimization of the modified objective functional (18). It turns out that depending on the initial guess for the parameters of the sequence, the optimal control algorithm reaches different local minima of (18) with distinct values. We have not studied this aspect in detail, but among the small sample of three different solutions we have obtained, the value of the minimum corresponding to the sequence shown in the top panel of Fig. 3 was the smallest. The application of each of the three sequences in MD simulations produced a long-lived stationary localized excitation like that shown in Fig. 3. From a practical viewpoint, the multiplicity of solutions suggests that it should be possible to impose additional constraints in the optimization algorithm to enhance the experimental feasibility of the resulting pulse sequence.

Although the results of Fig. 3 appear promising, they are sensitive to the system’s finite size as anticipated at the beginning of this section. For instance, applying the identical pulse sequence to a 42-particle system with FBCs yields no localization. In Fig. 3, the system’s ends are essential for the flow of vibrational energy towards the center to set up the target ILM. Hence this “direct” ILM excitation method is not suited for crystals, although it may have relevance for the excitation of anharmonic “local modes” in molecules such as benzene (C6H6). Indeed, ab initio MD simulations for benzene readily yield local modes, and the addition of an optimal control scheme may allow the prediction of optical waveforms for their creation.

#### 2 IR excitation

We have attempted to obtain similar results for the direct IR excitation of ILMs in finite systems, using a single chirped far-IR pulse of fixed width for the external force. Although we can reach a final state $`𝝃(t_f)`$ that is quite well localized, the localization does not persist for any significant time after the applied field ends. Evidently, the more restricted prescribed analytic time-dependence of the chirped Gaussian far-IR pulse reduces the ability for successful control.

### B “Indirect” excitation of ILMs in infinite lattices

The results of Sec. III A indicate that due to its size dependence, the direct excitation method is not a suitable approach for the creation of ILMs in crystals. However, our previous studies of the interrelation between ILMs and their associated ExM, as detailed in Refs.  and , suggest an alternative approach: because we have designed our model such that gap ILMs are related to the optically-active OZCM, we can create localized vibrations in the gap “indirectly” by optically driving the unstable OZCM. For this, our choice of a model in which the ILMs’ associated ExM is optically active is essential.

#### 1 ISRS excitation

To illustrate, we again use ISRS excitation. For the optimal control algorithm target state, we specify an OZCM of frequency $`\omega =0.98\omega _0`$, in an 8-particle lattice with StdPBCs. This particular system size is large enough to avoid computational complications due to second-neighbor interactions across the supercell boundaries, and it is sufficiently small to accelerate the optimal control scheme. The dynamics of the anharmonic OZCM are independent of the size of the StdPBC supercell, so that the resulting optimal fields apply to an infinite lattice. We choose this particular target OZCM on the basis of our RWA stability analysis of Sec. II C. At $`\omega =0.98\omega _0`$, the OZCM has an amplitude $`A=4.37\times 10^2R_0`$, a maximum instability growth rate $`\lambda _{\mathrm{max}}=2.78\times 10^2\omega `$, and the corresponding preferred perturbation wave vector is $`(k_p)_{\mathrm{max}}=0.624\pi /(2R_0)`$. Hence, this particular OZCM is expected to decay reasonably fast into ILM-like localized excitations with an easily discernible spacing of about $`6.4R_0`$. The advantage of choosing such an intermediate amplitude is clear from Fig. 2: at smaller amplitudes the expected spacing is larger, but the corresponding growth rate is smaller, and vice versa for larger amplitudes. The issue of the size of the growth rate is important, because we want to ensure that the decay of the unstable OZCM occurs on time scales where other processes, e.g., damping, which are not included in our description of the lattice dynamics, will not significantly alter the results. However, since an amplitude corresponding to $``$4% of the average equilibrium separation is quite large, we have repeated the optimization using as a target state an OZCM of smaller amplitude. The results for that case will be discussed at the end of this section.

Starting from rest at $`t=0`$, the system is driven with a sequence of 49 Gaussian pulses over a control interval of $`[0,50T_0]`$, as for the direct excitation detailed in Sec. III A 1. We again use weight factors $`\psi _n=16(m_n\omega _0/2)`$ and $`\psi _{n+N}=16/(m_n\omega _0/2)`$, but the simpler control task here allows us to increase the weight factor for the integrated square magnitude of the external force to $`\psi _{}=10`$, without affecting the ability to reach the target state. In contrast to the pulse sequence given in Fig. 3 for the finite chain, the optimal sequence for OZCM excitation is found to consist of pulses having nearly equal amplitudes. Accordingly, we simplified the control algorithm so as to vary the pulses’ common width, common amplitude, and individual pulse center times.

The top panel of Fig. 4 shows the resulting $`(t)`$, which consists of pulses of FWHM 32 fs and amplitude 0.013 eV/Å. This is an order of magnitude less than the largest amplitude for the direct ILM excitation in the finite 22-particle lattice of Fig. 3, and the equal amplitudes render this sequence qualitatively simpler. However, a closer look reveals important details, demonstrating how the control algorithm has globally optimized this pulse sequence. The solid line in the upper panel of Fig. 5 plots the position of a light particle in the OZCM as a function of time during the last fifth of the control interval. The thin vertical lines indicate the Gaussian pulse center times {$`t_i`$} of the external force. As expected for an efficient impulsive driving force, the pulse centers {$`t_i`$} coincide with the zero crossings of the particle position. In addition, we measured the “instantaneous” amplitude of the OZCM during the excitation by taking half the displacement difference between adjacent turning points of the motion of a light particle and assigning it the time of the intermediate zero crossing. Then we calculated the corresponding undriven RWA frequency. This is shown as a solid line in the lower panel of Fig. 5, while the circles indicate the frequencies $`2\pi /(t_{i+1}t_i)`$ corresponding to the spacings between adjacent pulses of the external force. The good agreement between these two quantities brings out the important (and somewhat hidden) aspect of the sequence of Fig. 4: the spacing between adjacent pulses varies through the sequence in such a way as to maintain resonant impulsive driving of the anharmonic OZCM, whose frequency decreases as its amplitude increases over the control interval. In Appendix A we show that the qualitative aspects of this behavior can be readily understood by doing the optimization for a sequence of $`\delta `$-function pulses.

Applying our Gaussian pulse sequence to a 40-particle StdPBC lattice in an MD simulation with the system initially at rest, we find that after the field ends at $`t_f=50T_0`$, the excited OZCM keeps vibrating with constant amplitude for about $`150T_0`$ until the perturbation due to accumulated computational round-off error triggers the instability and the OZCM decays into several localized excitations, as shown in the bottom panel of Fig. 4. The spatial array of localized excitations is not perfectly periodic because of the presence of instability perturbations of many wavevectors, with differing growth rates. Instead of relying on the clearly computer-dependent behavior of Fig. 4, we can provide the perturbation necessary to trigger the OZCM instability by including the effects of nonzero temperature. The bottom panel of Fig. 6 displays an MD simulation for the same pulse sequence, but with random initial velocities corresponding to a lattice temperature of 5 K. The presence of this perturbation triggers the OZCM decay much sooner. Of course the details of the MD results depend on the specific set of initial velocities, but for ten sets consistent with 5 K we find the same qualitative results as shown in Fig. 6: the ILM-like localized excitations resulting from the decay of the OZCM persist at fixed locations for several tens of vibrational periods and tend to move slowly through the lattice.

As mentioned earlier, we have repeated the indirect ILM excitation via ISRS with the same control period, the same number of pulses in the sequence, and identical weight factors, but using as a target state a less anharmonic OZCM at $`\omega =0.99\omega _0`$, with the corresponding amplitude $`A=3.03\times 10^2R_0`$. In this case, our optimal control algorithm yields a sequence of pulses with 32 fs FWHM and amplitude 0.009 eV/Å, compared with 32 fs and 0.013 eV/Å for the target OZCM at $`\omega =0.98\omega _0`$ and $`A=4.37\times 10^2R_0`$. The qualitative time evolution of the unstable OZCM excited using this optimal sequence in MD simulations with nonzero initial temperature is similar to that shown in Fig. 6, with the resulting ILM-like localized vibrations having smaller amplitudes. Moreover, as expected from the larger preferred instability length scale and smaller growth rate at this OZCM amplitude (see Fig. 2), the localized excitations are further apart from each other and it takes roughly 20$`T_0`$ longer than in Fig. 6 before a comparable degree of localization is reached.

#### 2 IR excitation

We have also studied the indirect ILM creation via OZCM excitation using IR, assuming for $`(t)`$ a single Gaussian pulse \[Eq. (30)\] of fixed width $`33T_0`$ (FWHM) and having a linearly chirped far-IR frequency over the control interval $`[0,100T_0]`$. Using the same OZCM target state and weight factors as detailed for the ISRS excitation above, the control scheme yields an optimal force with pulse amplitude 0.008 eV/Å and chirp rate $`2.0\times 10^7`$ fs-2. The corresponding $`(t)`$ is shown in the upper panel of Fig. 7, while the lower panel displays the results of applying this external force to a 40-particle StdPBC lattice in an MD simulation with initial random velocities corresponding to a lattice temperature of 5 K. Just as for the ISRS excitation of Fig. 6, the external force excites a slightly perturbed OZCM, which subsequently decays into ILM-like localized excitations. The above discussion concerning different sets of random velocities at the same temperature applies here as well.

Again, the optimization was repeated with the same control period, the same fixed pulse width, and identical weight factors, but using as a target state the smaller-amplitude OZCM at $`\omega =0.99\omega _0`$. In this case, our algorithm yields an optimal force with pulse amplitude 0.006 eV/Å and chirp rate $`4.2\times 10^8`$ fs-2. In analogy to the excitation of this OZCM with smaller amplitude via ISRS, the behavior in MD simulations with nonzero initial temperature is qualitatively similar to that for the larger-amplitude OZCM shown in Fig. 7, but exhibits the same differences as in the case of ISRS excitation: the resulting localized vibrations are spatially further apart and it takes about 20$`T_0`$ longer until a comparable degree of localization is reached.

## IV Discussion

### A Feasibility of the necessary external fields

Having theoretically demonstrated the creation of ILMs using “designer” external forces, we now discuss the experimental feasibility of the corresponding fields. First we consider the excitation via ISRS. Pulse shaping for ultrashort (13 fs) visible laser pulses has been demonstrated, with complicated final waveforms ranging from 12-pulse sequences with an overall Gaussian envelope and equal spacing to 6-pulse sequences with equal amplitudes and variable spacing. Furthermore, visible lasers producing ultrashort (18 fs) pulses with extremely large field magnitudes up to $``$270 V/Å are available. It is thus possible to produce the maximum field strengths of 1.22 V/Å and 0.38 V/Å necessary for the examples of Figs. 3 and 6, respectively, assuming a conservative value $`𝒫=2.5`$ Å2 for the polarizability derivative. However, the crucial experimental question is whether a given sample can tolerate such high electric fields in an experiment.

In a theoretical argument, Nelson and coworkers estimated the potential of single-pulse ISRS to excite large-amplitude anharmonic vibrations. Assuming pulses with 10 $`\mu `$J energy focused to 50 $`\mu `$m (FWHM) spot sizes, they predicted that a maximum phonon displacement of $`2\times 10^3`$ Å could be produced in the organic molecular crystal $`\alpha `$-perylene. For pulse widths of 70 fs, these parameters correspond to a field strength of about 0.5 V/Å. Hence, their estimation of the displacement produced by a single pulse and the field strength they considered reasonable are comparable with the corresponding quantities in our indirect ILM excitation via ISRS. Furthermore, they suggested that in some materials coherent vibrational displacements in the 0.1-1 Å range could possibly be achieved. These optimistic predictions were preceded by an actual experimental demonstration that single visible laser pulses of 70 fs width and 1 $`\mu `$J energy focused to 150 $`\mu `$m spot sizes impulsively excited coherent optic modes in $`\alpha `$-perylene. The field strengths, $`0.05V/\AA `$, used in this single pulse experiment were one order of magnitude below those assumed for the subsequent theoretical prediction. However, as discussed in a later experimental paper by the same group, it turned out not to be feasible to use the proposed larger field strengths, because they exceed the fairly low laser-induced breakdown threshold of $`\alpha `$-perylene. This sequence of publications highlights how a material’s laser-induced breakdown threshold can limit the potential of ISRS for the excitation of large-amplitude vibrations. However, it was pointed out in Ref.  that materials having a substantially larger laser-induced damage threshold than $`\alpha `$-perylene exist.

The question of breakdown thresholds in ultrashort pulse laser-solid interaction is not yet very well studied. The breakdown thresholds for alkali halides under irradiation by near-visible laser pulses with pulse widths down to 10 ps were experimentally determined to be about 0.2 V/Å (Ref. ), with a tendency for the thresholds to increase with increasing frequency and decreasing pulse width. Only few measurements for femtosecond pulses, such as the pulses used in our ISRS excitation studies, are available in the literature. For fused silica and the alkali fluorides, breakdown thresholds above 0.5 V/Å were obtained with 275 fs and 400 fs pulses for visible and near-visible frequencies. Values between 0.6 V/Å and 1.0 V/Å were measured for fused silica, sapphire, magnesium fluoride, and glass, using visible laser pulses having a width of 120 fs. Other experimental studies of fused silica obtained breakdown threshold fields of 3.8 V/Å, 3.0 V/Å, and 3.3 V/Å for visible pulses of 150 fs, 100 fs, and 55 fs width, respectively. We did not find measurements of the breakdown threshold in ZnS in the relevant short-pulse regime. However, an experiment on optical coatings made from ZnS (Ref. ) showed that its breakdown field strength for nanosecond pulses is comparable to that of magnesium fluoride, for which the short-pulse values are given above. From these experimental results for pulses that are still 2-20 times wider than the ones we used for the excitation of ILMs via ISRS in Secs. III A and III B, it appears that the maximum field strength for our direct ILM excitation of Fig. 3 for the 22-particle chain may be slightly too large to be realized in an experiment, while that for the indirect ILM excitation in the crystal lattice of Fig. 6 is below the values for breakdown. As discussed at the end of Sec. III B 1, the necessary external force magnitudes can be reduced by targeting an OZCM at a smaller amplitude, but this approach has limitations due to the competition between the time scales for the unstable OZCM decay and other processes in a crystal, e.g., damping. Moreover, the use of sequences with more pulses can also decrease the force magnitudes: according to the $`\delta `$-pulse approximation (see Appendix A), the force magnitude is inversely proportional to the number of pulses in the sequence. However, since a larger number of pulses requires longer control periods, the above caveat about the competition of time scales applies here as well.

Turning to IR excitation, we note that the maximum force amplitude in Fig. 7 corresponds to a field strength 0.008 V/Å, assuming that $`q=1.0e`$. Free-electron far-IR lasers produce picosecond pulses with intensities reported up to $`4\times 10^7`$ W/cm2 (Ref. ), corresponding to a field magnitude of 0.002 V/Å. Moreover, the frequency of free-electron far-IR lasers can be chirped at rates of $`9\times 10^9`$ fs-2 (Ref. ). These field magnitudes and chirp rates are within a factor of 5 and 20, respectively, of those used for the IR excitation of ILMs in Fig. 7. For this excitation mechanism, laser-induced breakdown should not play a role as a limiting factor, since the threshold of 0.2 V/Å obtained for alkali halides using near-IR pulses of 10 ps width  is well above the necessary maximum field magnitudes obtained here. As discussed at the end of Sec. III B 2, both the necessary field magnitudes and the chirp rates are reduced when a smaller amplitude OZCM is targeted.

We conclude that the fields necessary for the indirect excitation of ILMs via ISRS or far-IR as demonstrated in Sec. III B are reasonable and may be feasible in the near future. Among the approaches we have considered, excitation via ISRS seems more promising since lasers producing the necessary high field strengths are already available, although the experimental problem of laser-induced breakdown must be borne in mind.

Preliminary results indicate that the indirect excitation of ILMs can also be achieved when simpler analytic time dependences for the external force are used, without necessitating significantly larger force magnitudes. For the indirect ILM excitation via ISRS, we repeated the optimization procedure using pulse sequences constrained to have an equal variable spacing between the pulses as well as common variable widths and magnitudes and found that for the target OZCMs at $`\omega =0.98\omega _0`$ and $`\omega =0.99\omega _0`$ the required force magnitudes were 2.3% and 0.4% larger, respectively. Similarly, repeating the optimization for the indirect ILM excitation via IR using a single Gaussian pulse with unchirped variable frequency as well as fixed width and variable magnitude, we found the corresponding increases in the force magnitude to be 5.0% and 1.3%.

Another important experimental consideration is the robustness of the optimal fields. Considering the experimental limitations on the fidelity of shaped waveforms for ISRS, we note from Ref.  that pre-specified pulse amplitudes and widths were reproduced to within 10% and pulse positions to within 10 fs. Randomly perturbing the ISRS pulse parameters of the sequence shown in the top panel of Fig. 4 for the indirect ILM excitation in the lattice within these margins reveals that although the perturbed $`(t)`$ excites the OZCM to a slightly different amplitude in each of the ten cases considered, a decay into localized excitations always occurs. On the other hand, the ISRS pulse sequence shown in the top panel of Fig. 3 for the direct ILM excitation in the 22-particle chain is more sensitive to such infidelities. In each of ten cases of randomly perturbing the ISRS pulse parameters within the above margins, a localized excitation at the final time results, but only in one case does this excitation persist after the applied field is turned off. Decreasing the margins to a 5% error for pulse amplitudes and widths, and 5 fs for the pulse positions, the success rate increases to four out of ten.

### B Effects of nonzero initial temperature and damping

We have obtained our optimal external forces assuming that the system is initially at rest, but in an experiment thermal fluctuations will always be present and we should consider their effect. For the case of indirect ILM excitation, we have seen that the efficacy of the external field is enhanced by the presence of velocity perturbations due to an initial temperature of 5 K, since they serve to trigger the OZCM instability. At this low temperature, thermal fluctuations are small enough to act just as perturbations on the zero temperature dynamics. Accordingly, application of the zero-temperature optimal force initially excites a slightly perturbed OZCM which subsequently decays, as seen in the lower panels of Figs. 6 and 7, respectively.

If we increase the initial temperature, the situation changes. We have performed MD simulations with the ISRS pulse sequence of Fig. 4 for ten different sets of initial velocities appropriate to lattice temperatures of 77 K and 300 K. Already at 77 K the initial excitation can no longer be identified as a perturbed OZCM. Nevertheless, in all ten of the 77 K cases, the simulations result in well-localized ILM-like excitations, which appear sooner than for the 5 K case of Fig. 6. However, fewer localized excitations are observed for the same size lattice, and an increased background of thermally excited long-wavelength acoustic vibrations is present. This trend continues as the initial temperature is raised to 300 K. Very similar results are obtained when the optimal force for the far-IR excitation from Fig. 7 is used at higher initial temperatures. For comparison, we repeated the same nonzero initial temperature simulations, but with no external force, and we observed no significant localization of vibrational energy. Therefore, although the optimal external force obtained for zero initial temperature does not achieve its original goal of exciting an unstable OZCM when used at these elevated temperatures, it nevertheless produces localized vibrations.

Similar to the question of robustness with respect to infidelities in the pulse parameters, the optimal ISRS pulse sequence for the direct ILM excitation in the 22-particle chain is more sensitive to nonzero initial temperatures than is the sequence for indirect ILM excitation in the lattice. Performing MD simulations with the ISRS pulse sequence of Fig. 3 for ten different sets of initial velocities appropriate to lattice temperatures of 5 K, we find that a persisting localized excitation at the target site results in seven cases. Already at an initial temperature of 77 K the success rate drops to zero, although in some cases a localized excitation is created at a site other than the target site.

Vibrations in real crystals couple to other types of excitations and exhibit finite lifetimes – typically between 5 and 300 phonon periods for optical phonons. We have included this aspect by adding phenomenological velocity-dependent damping to our MD simulations. Anharmonicity, which contributes significantly to phonon lifetimes, is already treated explicitly in our calculations; thus we assume a small damping constant corresponding to an OZCM lifetime of $`100T_0`$. We repeated the optimization for the direct ILM excitation via ISRS of Fig. 3 after adding this damping, using the optimal pulse sequence with zero damping as initial guess. This results in a qualitatively very similar sequence of pulses, but with a larger maximum force magnitude of 0.24 eV/Å, compared with the earlier result 0.13 eV/Å for zero damping. Due to the presence of damping the amplitude of the resulting localized excitation decays away within a few tens of $`T_0`$ after the applied force ends. Similarly, when we include damping in the optimization of the indirect ILM excitation via ISRS of Fig. 6, we find that the common pulse amplitude increases from 0.013 eV/Å to 0.016 eV/Å. With $`T=0K`$ initial conditions, the amplitude of the resulting OZCM damps out before the accumulated computational round-off error can trigger the OZCM instability, but for an initial temperature of 5 K, we again observe a break-up of the OZCM into ILM-like localized excitations, whose amplitude subsequently damps out over few tens of OZCM periods. Additional calculations for the indirect ILM excitation via ISRS using a damping constant appropriate to a shorter OZCM lifetime of $`50T_0`$ yield a common pulse magnitude 0.019 eV/Å along with qualitatively similar MD results. Thus we conclude that it is possible to create ILMs with damping present, although the necessary force amplitudes are of course somewhat larger. Furthermore, detection of these ILMs would have to occur within few tens of OZCM periods after their creation.

### C Relevance for real crystals

While our 1D model incorporates some realistic features, such as standard interparticle potentials and the measured harmonic dispersion of ZnS, it is not a model of any real crystal. We have also considered 3D models of ZnS-structure crystals using standard two-body central potentials between atoms out to second neighbors. It is well known that the harmonic modes for $`𝐤`$ along $`<`$111$`>`$ map onto an effective 1D model involving the collective motion of (111) planes. We have shown that this mapping also occurs in the anharmonic case, for certain polarization directions. The resulting quasi-1D model, with effective anharmonic potentials between (111) planes undergoing collective motion with one transverse and one longitudinal degree of freedom, yields a representation of the actual motion in a 3D crystal and thus justifies our studies of purely 1D lattices. Applying again the ILM existence criterion of Ref. , we can choose the parameters of the quasi-1D model such that the ILMs are associated with the optically-active ExM. Hence, we adjusted our potential parameters so as to fit the measured harmonic phonon dispersion data of ZnS along $`<`$111$`>`$, as well as measured mode Grüneisen parameters. Within the RWA this model exhibits gap ILMs associated with the anharmonic version of the $`𝐤=0`$ TO phonon. These results suggests that to the extent that the central potentials used in this model capture the anharmonic properties of the real crystals, representative candidate materials for the indirect optical excitation of these ILMs would be ZnS, ZnSe, and the copper halides. We re-emphasize that these candidates for indirect optical creation of ILMs have two basic properties in common: (i) the $`𝐤=0`$ TO phonon in these materials is both first-order Raman and IR active, and (ii) the frequency of the $`𝐤=0`$ TO phonon is at the minimum of the TO branch along $`<`$111$`>`$, such that in conjunction with the fact that real potentials are dominated by soft anharmonicity, the criterion of Ref.  predicts the existence of gap ILMs associated with this $`𝐤=0`$ TO phonon.

## V Conclusion

In conclusion, our optical excitation studies demonstrate theoretically that suitably tailored laser radiation offers a promising route for the laboratory creation of vibrational ILMs. The time dependence of the fields is determined by an efficient optimal control algorithm, designed to produce waveforms consistent with the rapidly developing experimental capabilities in laser pulse shaping. The direct excitation of ILMs in finite systems was demonstrated for impulsive stimulated Raman scattering by a sequence of ultrashort laser pulses at THz repetition rates, with variable spacing between consecutive pulses. For periodic lattices, ILM creation was achieved indirectly via decay of the unstable associated ExM which is excited optically either via multiple-pulse ISRS as above or via a single far-IR pulse with a linearly chirped frequency. For the indirect excitation approach, it is essential to consider a lattice having ILMs that are associated with the optically-active ExM.

The direct ILM excitation via ISRS requires rather complex pulse sequences and laser field strengths that may just exceed the breakdown threshold of a given sample. Furthermore, this approach shows a dependence on the system size and is therefore not well suited for the creation of ILMs in crystal lattices, although it may be relevant for the excitation of local modes in molecules. Our studies show that the more advantageous means to excite ILMs in crystals is via the decay of their associated unstable anharmonic ExM, optically driven to a large amplitude. This approach not only requires smaller field strengths, but it succeeds even in the presence of thermal fluctuations and damping, and also with external forces constrained to have simpler analytic time dependences, i.e., pulse sequences with equal spacing between consecutive pulses for ISRS and a single pulse with constant frequency for IR excitation. Hence, our studies point to a potentially fruitful avenue for experimentally accessing the regime of large-amplitude anharmonic vibrational dynamics, which is very different than that for harmonic or weakly anharmonic systems.

###### Acknowledgements.

Supported by NSF Grant No. DMR-9510182. J. B. Page also gratefully acknowledges the Max Planck Institute for the Physics of Complex Systems, Dresden, Germany, for their hospitality and support during the completion of the manuscript.

## A Delta-approximation for ISRS excitation of the OZCM

In this Appendix we demonstrate that the main characteristics of the optimal pulse sequence for the ISRS excitation of the OZCM in our model lattice described in Sec. III B 1 can be reproduced using a simplification which was considered in Ref. , but in a different context. Here, the external force due to a sequence of short laser pulses is approximated by a sequence of $`L`$ delta functions

$$(t)=\underset{j=1}{\overset{L}{}}a_j\delta (tt_j)$$
(A1)
with individual positive amplitudes $`a_j`$ and pulse center times $`t_j`$. In addition, we know that the dynamics of the OZCM in a diatomic lattice map onto that of an effective anharmonic oscillator with mass $`m_{\mathrm{eff}}`$, displacement $`x`$ and momentum $`p`$. Including the external force (A1), the Hamiltonian for the driven case is given by

$$H_{\mathrm{eff}}=\frac{p^2}{2m_{\mathrm{eff}}}+V_{\mathrm{eff}}(x)x(t),$$
(A2)
where $`V_{\mathrm{eff}}(x)`$ is an effective anharmonic potential.

Starting with the effective oscillator initially at rest we want to excite it to a given final energy $`E_f`$ while keeping the necessary external force magnitude minimal. This requires an optimization of the parameters {$`a_j`$} and {$`t_j`$} for the external force (A1). A single delta pulse $`a_j\delta (tt_j)`$ boosts the momentum of the oscillator by $`a_j`$ at time $`t_j`$ and instantaneously changes the kinetic energy. Denoting the momentum of the oscillator immediately before and after the pulse by $`p_{j1}`$ and $`p_j`$, respectively, we can write the energy boost as

$$\mathrm{\Delta }E_j=\frac{\left(p_j^2p_{j1}^2\right)}{2m_{\mathrm{eff}}}=\frac{\left(a_j^2+2a_jp_{j1}\right)}{2m_{\mathrm{eff}}}.$$
(A3)
Evidently, for a given $`a_j`$ the largest possible energy transfer occurs when the pulse arrives exactly at the time when $`p_{j1}`$ is positive and maximal, i.e., at the zero crossings of the displacement of the oscillator where its momentum is positive. Knowing the effective potential, this simple result is sufficient to determine the optimal pulse center times for a given set of pulse amplitudes {$`a_j`$}. Furthermore, we now know that the energy of the oscillator after the complete pulse sequence is given by

$$E_L=\frac{1}{2m_{\mathrm{eff}}}\left(\underset{j=1}{\overset{L}{}}a_j\right)^2.$$
(A4)
We define an objective functional

$$J(\{a_j\})=\underset{j=1}{\overset{L}{}}a_j^2\lambda \left(E_LE_f\right),$$
(A5)
where the first term is just the integrated square magnitude of the external force (A1) and the second term constrains the energy $`E_L`$ after the pulse sequence to be equal to the desired final energy $`E_f`$ via introduction of a Lagrange parameter $`\lambda `$. Minimizing this objective functional with respect to the pulse amplitudes {$`a_j`$}, we find that the optimal pulse sequence consists of pulses with equal amplitudes

$$a_j=\frac{1}{L}\sqrt{2m_{\mathrm{eff}}E_f}.$$
(A6)
To summarize, the $`\delta `$-approximation predicts that the optimal pulse sequence for the OZCM excitation consists of pulses with equal amplitudes whose pulse center times coincide with the zero crossings of the motion of the atoms. This agrees with our numerical results from the full optimal control scheme, as discussed in Sec. III B 1 and illustrated in Fig. 5. However, the $`\delta `$-approximation underestimates the area under a single pulse by about 20%, when compared with the pulses in the numerically obtained optimal sequence. Hence, this simplified description correctly predicts the qualitative features of the optimal pulse sequence, but it cannot give quantitatively reliable results for the optimal pulse width and magnitude. These require knowledge of the actual motion of the particles over the duration of each pulse and are obtained by applying the full optimal control scheme.


format: text

text:
# 1 Introduction

## 1 Introduction

The $`SU(2)/U(1)`$ harmonic superspace (HSS) has been used to solve the $`D=4,N=2`$ off-shell constraints . The integrability interpretation of the dual version of the equations of motion has been also formulated for the $`N=2`$ supersymmetric Yang-Mills (SYM) theory in the framework of this superspace .

In the HSS-approach to the $`D=4,N=3`$ supersymmetry , the $`SU(3)/U(1)\times U(1)`$ harmonics have been used for the off-shell description of the corresponding SYM-theory . The $`N=3`$ SYM-equations in the ordinary superspace  have been transformed to the zero-curvature equations for the harmonic gauge connections, however, nobody tried earlier to solve these equations. Moreover, it has been shown that some set of analytic connections does not generate any solutions in the ordinary superspace .

We propose to solve first the harmonic zero-curvature equations in terms of the independent non-analytic superfield matrix $`v`$ ( bridge) describing the transform between different representations of the gauge group. Then we analyze the dynamical Grassmann (G-) analyticity equations for the composed on-shell harmonic connections constructed in terms of this superfield $`v`$.

It will be shown that the special non-covariant gauge choice for the matrix $`v`$ simplifies drastically the solution of all equations. The crucial feature of our gauge is the nilpotency of $`v`$ and the corresponding G-analytic harmonic connections. The Lorenz invariance of the SYM-equations is broken down to the $`SO(1,1)`$ subgroup in this gauge. We demonstrate that the 1-st order harmonic bridge equations with the nilpotent HSS-connections produce also the linear 2-nd order differential constraints for the bridge matrix.

Thus, we show that the $`N=3`$ SYM-equations in the ordinary superspace can be transformed to the linear solvable matrix differential equations in the harmonic superspace. This method allow us to analyze also solutions of the dimensionally reduced SYM-equations with 12 supercharges.

## 2 Harmonic representation of $`D=4,N=3`$ SYM constraints

We shall consider the non-covariant $`SO(1,1)\times SU(3)`$ representation of the $`D=4,N=3`$ superspace coordinates $`z`$

$`x^{1\dot{1}}x^{=|}=t+x^3,x^{2\dot{2}}x^==tx^3,x^{1\dot{2}}y=x^1+ix^2,x^{2\dot{1}}\overline{y}=x^1ix^2,`$
$`\theta _i^1=\theta _i^+,\theta _i^2=\theta _i^{},\overline{\theta }^{i\dot{1}}=\overline{\theta }^{i+},\overline{\theta }^{i\dot{2}}=\overline{\theta }^i.`$ (2.1)
where $`1,2,\dot{1},\dot{2}`$ are the $`SL(2,C)`$ indices, $`i=1,2`$ and 3 are indices of the fundamental representations of the group $`SU(3)`$. The $`SO(1,1)`$ weights of these coordinates are $`0,\pm 1,\pm 2`$, respectively.

Let us introduce the $`(4|6,6)`$-dimensional superspace gauge connections $`A(z)`$ and the corresponding covariant derivatives $``$

$`_\pm ^i=D_\pm ^i+A_\pm ^i,\overline{}_{i\pm }=\overline{D}_{i\pm }+\overline{A}_{i\pm },`$ (2.2)
$`_{=|}=_{=|}+A_{=|},_==_=+A_=,_y=_y+A_y,\overline{}_y=\overline{}_y+\overline{A}_y,`$
where the space-time derivatives and the Grassmann derivatives are considered.

The $`D=4,N=3`$ SYM-constraints  have the following reduced-symmetry form:

$`\{_+^k,_+^l\}=0,\{\overline{}_{k+},\overline{}_{l+}\}=0,\{_+^k,\overline{}_{l+}\}=2i\delta _l^k_{=|},`$ (2.3)
$`\{_+^k,_{}^l\}=\overline{W}^{kl},\{_+^k,\overline{}_l\}=2i\delta _l^k_y,`$ (2.4)
$`\{_{}^k,\overline{}_{l+}\}=2i\delta _l^k\overline{}_y,\{\overline{}_{k+},\overline{}_l\}=W_{kl},`$ (2.5)
$`\{_{}^k,_{}^l\}=0,\{\overline{}_k,\overline{}_l\}=0,\{_{}^k,\overline{}_l\}=2i\delta _l^k_=,`$ (2.6)
where $`W^{kl}`$ and $`\overline{W}^{kl}`$ are the gauge-covariant superfields constructed from the gauge connections. The equations of motion for the superfield strengthes follow from the Bianchi identities.

Let us analyze first Eqs.(2.3) which can be treated as integrability conditions for the positive-helicity connections. These conditions have the following pure gauge solution:

$$_+^k=g^1D_+^kg,\overline{}_{k+}=g^1\overline{D}_{k+}g,_{=|}=g^1_{=|}g.$$
(2.7)
Using the on-shell gauge condition $`g=1`$ we can obtain the simple general solution of Eqs.(2.3)

$$A_+^k=0,\overline{A}_{k+}=0,A_{=|}=0.$$
(2.8)
The analogous gauge conditions excluding the part of connections has been considered in Ref. for the self-dual $`4D`$ SYM-theory and in Ref. for the 10D SYM equations.

The $`SU(3)/U(1)\times U(1)`$ harmonics  parameterize the corresponding coset space. They form an $`SU(3)`$ matrix $`u_i^I`$ and are defined modulo $`U(1)\times U(1)`$ transformations

$$u_i^1=u_i^{(1,0)},u_i^2=u_i^{(1,1)},u_i^3=u_i^{(0,1)},$$
(2.9)
where $`i`$ is the index of the triplet representation of $`SU(3)`$. The complex conjugated harmonics have opposite $`U(1)`$ charges

$$u_1^i=u^{i(1,0)},u_2^i=u^{i(1,1)},u_3^i=u^{i(0,1)}.$$
(2.10)
These harmonics satisfy the following relations:

$`u_i^Iu_J^i=\delta _J^I,u_i^Iu_I^k=\delta _i^k,`$
$`\epsilon ^{ikl}u_i^1u_k^2u_l^3=1.`$ (2.11)
The $`SU(3)`$-invariant harmonic derivatives act on the harmonics

$`_J^Iu_i^K=\delta _J^Ku_i^I,_J^Iu_K^i=\delta _K^Iu_J^i,`$
$`[_J^I,_L^K]=\delta _J^K_L^I\delta _L^I_J^K.`$ (2.12)
We shall use the special $`SU(3)`$-covariant conjugation

$$\stackrel{~}{u_i^1}=u_3^i,\stackrel{~}{u_i^3}=u_i^1,\stackrel{~}{u_i^2}=u_i^2.$$
(2.13)
The corresponding conjugation of the harmonic derivatives is

$$\stackrel{~}{_3^1f}=_3^1\stackrel{~}{f},\stackrel{~}{_2^1f}=_3^2\stackrel{~}{f},$$
(2.14)
where $`f(u)`$ is a harmonic function.

We can define the real analytic harmonic superspace $`H(4,6|4,4)`$ with 6 coset harmonic dimensions and the following left and right coordinates:

$`\zeta =(\xi ^{=|},\xi ^=,\xi ,\overline{\xi },\theta _2^\pm ,\theta _3^\pm ,\overline{\theta }^{1\pm },\overline{\theta }^{2\pm }),`$
$`\xi ^{=|}=x^{=|}+i(\theta _3^+\overline{\theta }^{3+}\theta _1^+\overline{\theta }^{1+}),\xi ^==x^=+i(\theta _3^{}\overline{\theta }^3\theta _1^{}\overline{\theta }^1),`$
$`\xi =y+i(\theta _3^+\overline{\theta }^3\theta _1^+\overline{\theta }^1),\overline{\xi }=\overline{y}+i(\theta _3^{}\overline{\theta }^{3+}\theta _1^{}\overline{\theta }^{1+}),`$ (2.15)
where $`\theta _I^\pm =\theta _k^\pm u_I^k,\overline{\theta }^{\pm I}=\overline{\theta }^{\pm k}u_k^I`$.

The CR-structure in $`H(4,6|4,4)`$ involves the G-derivatives

$$D_\pm ^1,\overline{D}_{3\pm },$$
(2.16)
which commute with the harmonic derivatives $`D_2^1,D_3^2`$ and $`D_3^1`$.

These derivatives have the following explicit form in the analytic coordinates:

$`D_\alpha ^{(1,0)}D_\alpha ^1=_\alpha ^1/\theta _1^\alpha \overline{D}_\alpha ^{(0,1)}\overline{D}_{3\alpha }=_{3\alpha }/\overline{\theta }^{3\alpha },`$ (2.17)
$`D^{(2,1)}D_2^1=_2^1+{\displaystyle \frac{i}{2}}\theta _2^+\overline{\theta }^{1+}_{=|}+{\displaystyle \frac{i}{2}}\theta _2^+\overline{\theta }^1_\xi +{\displaystyle \frac{i}{2}}\theta _2^{}\overline{\theta }^{1+}\overline{}_\xi +{\displaystyle \frac{i}{2}}\theta _2^{}\overline{\theta }^1_=`$
$`\theta _2^+_+^1\theta _2^{}_{}^1+\overline{\theta }^{1+}\overline{}_{2+}+\overline{\theta }^1\overline{}_2,`$
$`D^{(1,2)}D_3^2=_3^2+{\displaystyle \frac{i}{2}}\theta _3^+\overline{\theta }^{2+}_{=|}+{\displaystyle \frac{i}{2}}\theta _3^+\overline{\theta }^2_\xi +{\displaystyle \frac{i}{2}}\theta _3^{}\overline{\theta }^{2+}\overline{}_\xi +{\displaystyle \frac{i}{2}}\theta _3^{}\overline{\theta }^2_=`$
$`\theta _3^+_+^2\theta _3^{}_{}^2+\overline{\theta }^{2+}\overline{}_{3+}+\overline{\theta }^2\overline{}_3,`$ (2.18)
$`D^{(1,1)}D_3^1=_3^1+i\theta _3^+\overline{\theta }^{1+}_{=|}+i\theta _3^+\overline{\theta }^1_\xi +i\theta _3^{}\overline{\theta }^{1+}\overline{}_\xi +i\theta _3^{}\overline{\theta }^1_=`$
$`\theta _3^+_+^1\theta _3^{}_{}^1+\overline{\theta }^{1+}\overline{}_{3+}+\overline{\theta }^1\overline{}_3,`$ (2.19)
where $`_{=|}=/\xi ^{=|},_==/\xi ^=,_\xi =/\xi `$ and $`\overline{}_\xi =/\overline{\xi }`$.

It is crucial that we start from the specific gauge conditions (2.8) for the $`N=3`$ SYM-connections which break $`SL(2,C)`$ but preserve the $`SU(3)`$-invariance. Consider the harmonic transform of the covariant Grassmann derivatives in this gauge using the projections on the $`SU(3)`$-harmonics

$`_+^Iu_i^ID_+^i=D_+^I,\overline{}_{I+}u_I^i\overline{D}_{i+}=\overline{D}_{I+},\{D_+^I,\overline{D}_{K+}\}=2i\delta _K^I_{=|},`$ (2.20)
$`_{}^Iu_i^I_{}^i=D_{}^I+𝒜_{}^I,\overline{}_Iu_I^i_i=\overline{D}_I+\overline{𝒜}_I,`$ (2.21)
where the harmonized Grassmann connections $`𝒜_{}^I`$ and $`\overline{𝒜}_I`$ are defined.

The $`SU(3)`$-harmonic projections of the superfield constraints (2.4-2.6) can be derived from the basic set of the $`N=3`$ super-integrability conditions for two components of the harmonized connection:

$`D_+^1𝒜_{}^1=\overline{D}_{3+}𝒜_{}^1=D_+^1\overline{𝒜}_3=\overline{D}_{3+}\overline{𝒜}_3=0`$ (2.22)
$`D_{}^1𝒜_{}^1+(𝒜_{}^1)^2=0,\overline{D}_3\overline{𝒜}_3+(\overline{𝒜}_3)^2=0,`$
$`D_{}^1\overline{𝒜}_3+\overline{D}_3𝒜_{}^1+\{𝒜_{}^1,\overline{𝒜}_3\}=0.`$ (2.23)
All projections of the SYM-equations can be obtained by the action of the harmonic $`SU(3)`$ derivatives $`D_K^I`$ on these basic conditions.

These Grassmann zero-curvature equations have the very simple general solution

$$𝒜_{}^1(v)=e^vD_{}^1e^v,\overline{𝒜}_3(v)=e^v\overline{D}_3e^v,$$
(2.24)
where the bridge $`v`$ is the gauge-Lie-algebra valued (5,5)-superfield matrix satisfying the additional constraint

$$(D_+^1,\overline{D}_{3+})v=0.$$
(2.25)
Consider the gauge transformations of the bridge

$$e^ve^\lambda e^ve^{\tau _r},$$
(2.26)
where $`\lambda H(4,6|4,4)`$ is the (4,4)-analytic matrix parameter, and the parameter $`\tau _r`$ does not depend on harmonics. The matrix $`e^v`$ realizes the harmonic transform of the gauge superfields $`A_\pm ^k,\overline{A}_{k\pm }`$ in the central basis (CB) to the harmonic gauge superfields of the analytic basis (AB) using the analytic $`\lambda `$-representation of the gauge group.

The bridge $`v`$ defines the on-shell harmonic connections of the analytic basic

$$V_K^I(v)=e^vD_K^Ie^v$$
(2.27)
which satisfy by construction the harmonic zero-curvature equations, for instance,

$`D_2^1V_3^2D_3^2V_2^1+[V_2^1,V_3^2]V_3^1=0,`$
$`D_2^1V_3^1D_3^1V_2^1+[V_2^1,V_3^1]=0.`$ (2.28)
The dynamical SYM-equations in the bridge representation (2.24) are reduced to the harmonic differential conditions for the basic Grassmann connections:

$$(D_2^1,D_3^2,D_3^1)(𝒜_{}^1(v),\overline{𝒜}_3(v))=0,$$
(2.29)
which are equivalent to the following set of the dynamic G-analyticity equations:

$$(D_{}^1,\overline{D}_3)(V_2^1(v),V_3^2(v),V_3^1(v))=0.$$
(2.30)
The positive-helicity analyticity conditions

$$(D_+^1,\overline{D}_{3+})(V_2^1(v),V_3^2(v),V_3^1(v))=0$$
(2.31)
are satisfied automatically for the bridge in the gauge (2.25). Stress that the analyticity equations in the bridge representation describe all SYM-solutions without a loss of generality.

The inverse harmonic transform determines the harmonic-independent on-shell gauge solutions in the ordinary superspace

$`𝒜_{}^1(v)A_{}^i(v)=(u_1^i+u_2^iD_1^2+u_3^iD_1^3)e^vD_{}^1e^v`$
$`\overline{𝒜}_3(v)\overline{A}_i(v)=(u_i^3u_i^1D_1^3u_i^2D_2^3)e^v\overline{D}_3e^v.`$ (2.32)
By construction, these superfields satisfy the $`D=4,N=3`$ CB-constraints (2.4-2.6) and the harmonic differential conditions

$$D_K^I(A_{}^i(v),\overline{A}_i(v))=0,$$
(2.33)
if the relations (2.29) are fulfilled.

It is useful to calculate the Grassmann connection $`𝒜_{}^1(v)`$ in terms of the basic analytic matrices

$$e^vD_{}^1e^v=b^1\theta _1^{}(b^1)^2+\overline{\theta }^3\left(d_3^1\frac{1}{2}\{b^1,c_3\}\right)+\theta _1^{}\overline{\theta }^3\left([b^1,d_3^1]+\frac{1}{2}[c_3,(b^1)^2]\right).$$
(2.34)
The conjugated connection $`\overline{𝒜}_3(v)`$ can be constructed analogously.

## 3 Harmonic-superspace equations of motion

Using the $`\lambda `$-transformations (2.26) one can gauge away the first terms in the Grassmann decomposition of the bridge $`v`$ and choose the non-supersymmetric nilpotent gauge for this (5,5)-superfield:

$`v=\theta _1^{}b^1+\overline{\theta }^3c_3+\theta _1^{}\overline{\theta }^3d_3^1,v^2=\theta _1^{}\overline{\theta }^3[c_3,b^1],`$ (3.1)
$`e^v=Iv+{\displaystyle \frac{1}{2}}v^2=I\theta _1^{}b^1\overline{\theta }^3c_3+\theta _1^{}\overline{\theta }^3({\displaystyle \frac{1}{2}}[c_3,b^1]d_3^1),`$ (3.2)
where the fermionic matrices $`b^1,c_3`$ and the bosonic matrix $`d_3^1`$ are (4,4)-analytic functions of the coordinates $`\zeta `$ (2.15).

The restrictions on the bridge matrix in the gauge group $`SU(n)`$

$$\text{Tr}v=0,v^{}=v$$
(3.3)
are equivalent to the following relations for the analytic matrices $`b^1,c_3`$ and $`d_3^1`$

$$\text{Tr}b^1=0,\text{Tr}c_3=0,\text{Tr}d_3^1=0,(b^1)^{}=c_3,(d_3^1)^{}=d_3^1.$$
(3.4)
It is useful to construct the on-shell superfield strength in the analytic basis

$$\overline{W}^{12}(b^1)=D_+^2b^1+\theta _2^{}D_{}^2D_+^2b^1+\theta _2^{}[D_1^2b^1,D_+^2b^1]+\theta _1^{}\theta _2^{}[D_{}^2b^1,D_+^2b^1],$$
(3.5)
which satisfies the non-Abelian conditions of Grassmann and harmonic analyticities.

The dynamical G-analyticity equation (2.30) in the gauge (3.1) is equivalent to the following harmonic differential bridge equation:

$`V_2^1(v)=\theta _2^{}b^1v_2^1,(v_2^1)^2=0`$ (3.6)
$`D_2^1e^v=e^vv_2^1.`$ (3.7)
where the analytic nilpotent representation for the on-shell harmonic connection arises.

The 2-nd on-shell harmonic connection is also nilpotent

$`V_3^2(v)=\overline{\theta }^2c_3v_3^2,`$ (3.8)
$`D_3^2e^v=e^vv_3^2.`$ (3.9)
Underline that the nilpotency of the analytic parts of the harmonic connections $`V_2^1`$ and $`V_3^2`$ (2.27) follows directly from the gauge condition (2.31), and the analyticity of these connections requires the additional harmonic restrictions on the functions $`b^1,c_3`$ and $`d_3^1`$.

The 1-st bridge equation (3.7) generates the following nonlinear equations for the (4,4)-analytic matrices:

$`D_2^1b^1=\theta _2^{}(b^1)^2,`$ (3.10)
$`D_2^1c_3=\theta _2^{}(d_3^1+{\displaystyle \frac{1}{2}}\{b^1,c_3\}),`$ (3.11)
$`D_2^1d_3^1={\displaystyle \frac{1}{2}}\theta _2^{}\left([d_3^1,b^1]+{\displaystyle \frac{1}{2}}[(b^1)^2,c_3]\right).`$ (3.12)
It is important that all these equations contain the nilpotent element $`\theta _2^{}`$ in the nonlinear parts, so they can be reduced to the set of linear iterative equations.

The 2-nd bridge equation (3.9) gives us

$`D_3^2b^1=\overline{\theta }^2(d_3^1+{\displaystyle \frac{1}{2}}\{b^1,c_3\}),`$ (3.13)
$`D_3^2c_3=\overline{\theta }^2(c_3)^2`$ (3.14)
$`D_3^2d_3^1={\displaystyle \frac{1}{2}}\overline{\theta }^2\left([d_3^1,c_3]+{\displaystyle \frac{1}{2}}[b^1,(c_3)^2]\right).`$ (3.15)
It is useful to derive the following relations:

$`\theta _2^{}D_2^1(b^1,c_3,d_3^1)=0,\overline{\theta }^2D_3^2(b^1,c_3,d_3^1)=0,`$ (3.16)
$`\theta _2^{}D_3^2b^1+\overline{\theta }^2D_2^1c_3=\theta _2^{}\overline{\theta }^2\{b^1,c_3\},`$ (3.17)
$`\overline{\theta }^2D_2^1c_3\theta _2^{}D_3^2b^1=2\theta _2^{}\overline{\theta }^2d_3^1,`$ (3.18)
$`\theta _2^{}D_3^1b^1+\overline{\theta }^1D_2^1c_3=\theta _2^{}\theta _3^{}(b^1)^2+\theta _2^{}\overline{\theta }^1\{b^1,c_3\}.`$ (3.19)
The bridge equations with the nilpotent analytic connections yield the following 2-nd order linear conditions for the coefficient functions

$$D_2^1D_2^1(b^1,c_3,d_3^1)=0,D_3^2D_3^2(b^1,c_3,d_3^1)=0.$$
(3.20)
Now we consider the iterative procedure of solving the basic non-Abelian harmonic differential equations for the (4,4) analytic matrices $`b^1`$ and $`c_3`$ using the partial decomposition in the Grassmann variables $`\theta _2^{},\theta _3^{},\overline{\theta }^1,\overline{\theta }^2`$. The matrix (4,0) coefficients of this decomposition have dimensions $`1/2l5/2`$

$`b^1=\beta ^1+\theta _2^{}B^{12}+\theta _3^{}B^{13}+\overline{\theta }^1B^0+\overline{\theta }^2B_2^1+\theta _2^{}\theta _3^{}\beta ^0+\theta _2^{}\overline{\theta }^1\beta ^2+\theta _3^{}\overline{\theta }^2\beta _2^{13}`$
$`+\theta _3^{}\overline{\theta }^1\beta ^3+\overline{\theta }^1\overline{\theta }^2\beta _2+\theta _2^{}\overline{\theta }^2\eta ^1+\theta _2^{}\theta _3^{}\overline{\theta }^1B^{23}+\theta _3^{}\overline{\theta }^1\overline{\theta }^2B_2^3+\theta _2^{}\theta _3^{}\overline{\theta }^2C^{13}`$
$`+\theta _2^{}\overline{\theta }^1\overline{\theta }^2C^0+\theta _2^{}\theta _3^{}\overline{\theta }^1\overline{\theta }^2\eta ^3.`$ (3.21)
The first iterative (4,0) equations $`(l=1/2)`$ are linear and homogeneous

$$\widehat{D}_2^1\beta ^1=\widehat{D}_3^2\beta ^1=0.$$
(3.22)
The next harmonic iterative equations for the (4,0) components with $`l1/2`$ can be resolved on the each step via functions of the highest dimensions or their derivatives. These linear equations contain nonlinear sources constructed from the solutions of the previous equations. Note that some (4,0) iterative equations are pure algebraic relations which reduce the number of independent functions, for instance,

$`B_2^1=\widehat{D}_2^1B^0{\displaystyle \frac{i}{2}}\theta _2^+_\xi \beta ^1,B^{12}=\widehat{D}_3^2B^{13}+{\displaystyle \frac{i}{2}}\overline{\theta }^{2+}\overline{}_\xi \beta ^1,`$ (3.23)
$`\beta ^2=\widehat{D}_3^2\beta ^3{\displaystyle \frac{i}{2}}\overline{\theta }^{2+}\overline{}_\xi B^0,\beta _2^{13}=\widehat{D}_2^1\beta ^3+{\displaystyle \frac{i}{2}}\theta _2^+_\xi B^{13}.`$ (3.24)
The finite set of the linear harmonic differential equations for the independent (4,0) functions $`\beta ^1,B^0,B^{13},\beta ^3\mathrm{}`$ can be, in principle, explicitly solved. Thus, the $`SU(3)/U(1)\times U(1)`$ harmonic method together with the simple gauge conditions allow us to transform the $`N=3`$ superfield SYM-constraints to the harmonic differential equations which are equivalent to the finite set of the iterative solvable linear equations.

The authors are grateful to E.A. Ivanov for the discussions.

This work is supported by the Votruba-Blokhintsev programme in Joint Institute for Nuclear Research. The work of B.Z. is partially supported also by the grants RFBR-99-02-18417, RFBR-DFG-99-02-04022 and NATO-PST.CLG-974874.