text:
# Multiparticle Biased DLA with surface diffusion: a comprehensive model of electrodeposition

## I Introduction

Quasi-two-dimensional (quasi-2D) electrochemical deposition (ECD)  has become one of the most widely studied pattern forming processes since its recognition as a paradigm of non-local, non-equilibrium growth processes . Within this general context, a great deal of work has been devoted in the past fifteen years to experimental and theoretical studies of quasi-2D ECD. A first group of works deals mainly with pattern formation, its main results concerning “phase diagrams” of morphologies , ECD as Laplacian growth process , dynamic morphological transitions , etc. All these studies aim to understanding the principles underlying the rich variety of morphologies observed, ranging from dendritic to fractal. In addition to this line of research, there is a second one  whose main interest is the existence of universality and scale invariance in the roughness of the deposits produced . From all this and related research, it is now believed that complex structures with different morphologies arise from quasi-2D ECD due to the interplay of different transport mechanisms, such as cation diffusion, electromigration, fluid convection, and surface diffusion (SD) . However, the combined effect of all these factors leads to a very complex process, and it is becoming increasingly apparent that ECD is not well understood yet. In particular, the detailed role of surface diffusion (SD) is still an open question that hinders our understanding of both the morphologies and the scaling of ECD aggregates.

Much of the work mentioned in the above paragraph has been motivated by the quest to find a universal model to help understand ECD phenomena. The first model formulated with that purpose was the famous computer algorithm known as Diffusion Limited Aggregation (DLA) , in which a particle diffuses on a lattice and attaches to the growing aggregate at the place where it first hits. It is not difficult to observe (see  and references therein) that this simple model represents the zero concentration, quasi-static limit of ECD. Therefore, its validity as a general description of ECD is rather restricted because it does not include most of the effects involved in the process. However, DLA has played a seminal rôle as a source of inspiration both for continuum approaches  —which predict some high-current properties but take into account neither the influence of the applied voltage nor the electrolyte concentration— and for more sophisticated computer models, basically modifications of DLA (see, e.g.,  and also the paragraph below), which are more or less phenomenological and concentrate on changes in morphology, thus being unable to explain the underlying mechanisms yielding those patterns.

In this paper, we report on the results of detailed numerical studies of Multiparticle Biased Diffusion Limited Aggregation (MBDLA) supplemented with SD. MBDLA is a model in the family of multiparticle DLA models , in which a finite number of random walkers, possibly with constant concentration, is introduced instead of the single walker of DLA. Thus, the excluded volume interaction among the walkers leads to several of the effects neglected in DLA. As its main ingredient, MBDLA includes, in addition, a preferential bias (which had been first studied in the context of single-particle, DLA models by Meakin ) of the walkers towards the cathode to mimic the electric field: In this form, the model was successfully introduced in  to study the influence of the applied electric field on the composition of magnetic, amorphous CoP alloys grown by ECD at constant current. The main virtue of MBDLA is that it is a mesoscopic model embedded on a two-dimensional square lattice, but it reproduces the mean fraction of Co and P in two-cation species ECD, as well as the qualitative morphology of the product electrodeposits. In fact, the agreement between MBDLA and ECD experiments  is quantitative, as the electrical current intensity and the experiment time can be directly related to the simulation parameters . Therefore, we are confident that MBDLA is a good starting point to study the relevance of SD in ECD and, specifically, its influence on the shape of the aggregates and their dynamic scaling. Scaling properties of MBDLA without SD were briefly reported in .

The report of our results is organized as follows: We describe our model in Sec. II, where a brief introduction to the physics and chemistry of ECD is followed by a detailed account of the rules governing MBDLA. Section III reports on our numerical results, such as morphological patterns and roughness scaling. After physically showing that SD has to be included, we introduce three different rules for SD are carefully considered and compared to experiments, allowing to identify the proper way to introduce SD in the model. Finally, we conclude in Sec. IV with a discussion of our results which will allow us to suggest a reasonably approximate picture of ECD phenomena. A few technical details about one the rules for SD are given in an Appendix.

## II The Model

### A Basic facts about ECD

Prior to describing in detail what MBDLA is, and in order to motivate and to better understand the model rules, we will briefly summarize the basic physics and chemistry of ECD, by collecting the equations commonly accepted to govern its main features (see, e.g.,  for further details). Generally speaking, ECD experiments involve two species, named cations and anions, moving in an incompressible viscous fluid. In very many cases, ECD takes place in quasi-2D cells with parallel electrodes. The cations move towards the cathode and the anions towards the anode. The basic equations for the concentrations of both species are as follows:

$`{\displaystyle \frac{C}{t}}`$ $`=`$ $`J_c,`$ (2)
$`{\displaystyle \frac{A}{t}}`$ $`=`$ $`J_a,`$ (3)
$`𝐉_c`$ $`=`$ $`D_cC+\mu _c𝐄C+𝐯C,`$ (4)
$`𝐉_a`$ $`=`$ $`D_aA\mu _a𝐄A+𝐯A,`$ (5)
where $`C`$ and $`A`$ are the cation and anion concentration respectively, $`D_{c,a}`$ the cationic and anionic diffusion coefficients, $`\mu _{c,a}`$ their mobilities, $`𝐯`$ the fluid velocity field, and $`𝐄`$ the electric field along the cell. The latter is related to cation and anion concentration via the Poisson equation

$$E=^2\varphi =e(z_cCz_aA)/\epsilon ,$$
(6)
where $`\varphi `$ is the applied potencial, $`ez_c`$ and $`ez_a`$ are the cation and anion electric charges, respectively, and $`\epsilon `$ the dielectric permittivity of the fluid. Generally speaking, matter balance across the interface leads to an interface velocity proportional to the flux of cations and therefore, in the absence of any other limiting process, proportional to the current density as well. In addition, except for the region close to the cathode, we may assume electro-neutrality , which in turn implies that the cation mean velocity is constant. The incompressible Navier-Stokes equations determine the velocity $`𝐯`$ of the solvent. Fluid convection is always present in ECD experiments, but in many instances it can be small enough to be safely neglected, as has been shown in actual experiments .

When the cations arrive at the cathode, they reduce irreversibly and an aggregate of neutral particles begins to grow. The particles on the surface aggregate are transported all along it due to local chemical potential gradients. The resulting particle current conserves the number of particles on the surface and is given approximately by by (see  for a detailed discussion):

$$𝐉_s_s\kappa ,$$
(7)
where $`𝐉_s`$ is the particle current along the surface, $`\kappa `$ the interface local curvature at each site, and $`_s`$ the gradient taken along the surface. Roughly speaking, SD tends to reduce the interface local curvature. Finally, we note that the mean concentration of charge carriers in the bath is constant as new cations are formed at the anode upon arrival of the anions .

### B Definition and rules of MBDLA

In this section we will define MBDLA through its evolution rules, for which we take into account the physical equations presented in the previous section. At this point, we do not consider SD, whose need will be justified in the next Section, and consequently we postpone the discussion of the rules to implement SD as well. Thus, MBDLA is a cellular automaton defined on a two-dimensional square lattice of horizontal dimension $`L_x`$ and vertical dimension $`L_y`$ (with lateral periodic boundary conditions and reflective boundary condition at the top; for the conditions at the bottom, see below), in which a number of random walkers (cations) are randomly distributed with concentration $`c`$. The bottom of the lattice is chosen to be the cathode. We do not consider the anion dynamics, but we implicitly introduce it by the creation of particles and by charge electro-neutrality .

The initial condition evolves in time as follows: Every time step a walker is chosen and moved to one of its four neighboring sites with probabilities taken from a finite differences scheme of Eqs. (2) and (4) : probability $`1/(4+p)`$ to move either left, right or upwards, and probability $`(1+p)/(4+p)`$ to move down, i.e., towards the cathode. The parameter $`p`$ is referred to as the bias; in galvanostatic conditions it can be quantitatively related to the electric current density in the physical system as shown in . Let us stress here that our present choice for the probabilities is different from that reported in  and , but we have checked that the results hardly differ with those presented in this paper. The main reason for this new election is that, with the new rules, the bias $`p`$ ranges from $`0`$ to $`\mathrm{}`$, that is, from pure multiparticle DLA to ballistic deposition, whereas the rules in the mentioned references allow for a range in $`p`$ from $`0`$ to $`0.25`$, and the ballistic deposition limit cannot be reached (although $`p=0.25`$ is rather close already, see ). After a destination site has been chosen, the particle moves if that node of the lattice is empty; if not, we select another particle and repeat the destination selection procedure. Once the particle has been moved, if the new position has any nearest neighbor site belonging to the aggregate, the present position of the walker is added to the aggregate (the cathode or bottom boundary at the initial stage) with probability $`s`$ (and is able to diffuse over the aggregate surface if that aggregate position has just one nearest neighbor belonging to it, see the following section); otherwise it stays there (and is able to move again) with probability $`1s`$. We term $`s`$ the sticking probability; it is related to the chemical activation energy the cation needs to stick to the aggregate. As particles are added to the aggregate, other ones are created at the top of the lattice keeping the mean cation concentration $`c`$ constant, which in fact simulates an infinitely high system (experimentally this means that the distance between electrodes is much larger than their lateral dimension); consequently, the flux of particles is constant at every stage of the simulation.

As we have already pointed out, the model parameters are related to the physical factors influencing the problem. Indeed, the choice of jump probabilities for the random walkers in the bath provides a recurrence relation which is a discretized version of the continuous equations (2) and (4). Therefore, the drift velocity $`\mu _𝐜𝐄`$ is proportional to the bias $`p`$. When a finite number of walkers is considered with concentration $`c`$, we must take into account the excluded volume, so the effective diffusion coefficient and the effective drift velocity in the simulations are proportional to $`1c`$ (in a mean field approach) . It is important to note that when $`c0`$, i.e., the bath is formed by one particle alone (as in DLA), the aggregate develops tall branches which grow at the expense of short ones due to screening effects. Therefore, in the low current limit a morphological instability appears that is not always present in ECD experiments. The finite concentration and the hard core interaction among random walkers simulates the cation pressure on the aggregate, so $`c`$ is an essential ingredient in the understanding of the formation of electrodeposits and to prevent these instabilities (of Laplacian character) from dominating the whole growth process.

One important task is the definition of simulation time step. In , comparison with the experiments in  allowed to show that the physical time and the simulation time measured in number of Monte Carlos trials were simply proportional to each other. For this reason, we have stuck to the definition of the time step in  as a Monte Carlo trial, i.e., the time needed for a particle to jump, both if the particle does jump or if it does not. Notwithstanding, we have tried other time steps definitions, such as the Monte Carlo step being defined as the mean time for every random walker to jump at least once, but the results are basically the same. Some authors define the time step for solid-on-solid growth models as the mean time needed to complete an aggregate layer, but as we will show below, ECD electrodeposits do not grow with constant velocity, and therefore the mean interface height does not grow linearly with time. We thus believe that, in the ECD context, this time unit would be rather artificial and hence we have not used it. In fact, as we will show below, the work reported in this paper provides further evidence in favor of our choice (see the discussion of the experiments in  in Sec. III B below).

## III Numerical Results

### A Morphologies

We begin the summary of our results by discussing the morphologies generated by MBDLA with and without SD and the influence of the different rules for SD on them. In addition, we want to compare our computer generated morphologies to the available experimental data. We take as a reference the comprehensive experimental work of Trigueros et al. , who reported a systematic experimental study of different growth regimes at constant applied voltage conditions. Their work gave rise to a diagram of morphologies divided into different regions in which similar morphologies were obtained as a function of the applied voltage and the electrolyte concentration. It is important to realize that, in galvanostatic conditions, there is no linear correspondence between voltage and electric current of ions, and therefore, comparison between our morphologies and those reported by these authors can only be qualitative. No similar taxonomy work has been performed for constant current conditions. Although the diagram in  is quite complex, it encloses a full variety of morphologies under the label compact. Some authors  have studied electrodeposits systems within this regime, and hereafter we will also refer to them. Finally, a recent work by Schilardi et al.  provides exhaustive information on the asymptotic ECD regimes, which have not been considered anywhere else; hence, their research will also be compared to ours along the paper.

As we have already mentioned, from the model perspective we can compare the bias $`p`$ with the electric current density , and $`c`$ with the electrolyte concentration, even though the two latter magnitudes are not exactly coincident, i.e., an electrolyte concentration equal to $`0.1`$ M does not mean $`c=0.1`$. We will see below that the results are not very sensitive to the specific value of $`c`$ insofar it is not very small, and thus the difference between actual and model concentrations is not very relevant. The sticking probability, $`s`$, and the diffusion parameters, namely $`l`$, $`\lambda `$ and $`r`$ (or equivalently $`\tau _d`$), cannot be directly tuned in an experiment, although it is reasonable to expect that changes in the experimental conditions will correspondingly modify these parameters. How much are they modified is something we will learn through our computer simulations.

#### 1 Bias vs. Sticking Probability without SD

Figure 1 shows a diagram of morphologies obtained with $`0p5`$ and $`0.01s0.5`$ without SD, with a particle concentration $`c=0.05`$. We have included these results for two reasons: First, there has been no previous report on MBDLA morphologies, except for a brief discussion in ; and, second, we need to discuss them in order to understand later what is the effect of SD on MBDLA morphologies. It is clear from Fig. 1 that increasing the bias or decreasing the sticking probability yields denser aggregates, the ones obtained for $`p=0`$ and $`s=1`$ (bottom right) being multiparticle DLA-like as expected (compare to ). This phenomenon is related to the stabilizing effect of the parameters $`p`$ and $`s`$, which can be theoretically demonstrated . Indeed, the higher the values of $`p`$, the larger the flux of particles reaching the interface in the direction perpendicular to the cathode. This reduces the probability for a cation to stick laterally to a branch and the screening effects due to the Laplacian field. On the other hand, the electric field combined with the reduction of the sticking probability tends to fill the interface valleys. This first result, namely the fact that increasing the electric current leads to denser aggregates, is similar to the results reported by Trigueros et al. , who observed densification of the aggregates with increasing applied voltage. In particular, we can qualitatively compare the morphological changes obtained by varying the bias $`p`$ for a fixed $`s=0.5`$ in Fig. 1, with those provided by experimental voltage variations (see Fig. 2 in ). We conclude that high voltages (or in general, high density currents) yield denser aggregates. So the bias $`p`$ is an essential ingredient in any realistic ECD model.

As a second step in our study, we have monitored other relevant quantities which in turn can be experimentally measured, in order to obtain additional information aside from qualitative morphological comparisons. Figure 2 shows the local concentration of particles in the bath, still without SD, at equal time intervals. We have plotted the concentration profiles in the stationary regime, i.e., after the instability occurs (see below). Thus, the mean number of attached particles per unit time (or equivalently, the mean interface velocity) is constant. Léger et al.  have reported experimental evidence consistent with this stationary behavior (see e.g. Fig. 5 in ). We thus see that MBDLA agrees well with their findings, i.e., the stationary concentration of particles in the bulk obeys approximately the equation :

$$C(z)=c_a+(c_0c_a)e^{(zz_0)u/D_c},$$
(8)
where $`z`$ is the vertical coordinate, $`z_0`$ is the interface mean position, $`c_a`$ is the concentration at the anode, $`c_0`$ the concentration at the interface, $`u=\mu _c𝐄`$, and $`D_c`$ is the bulk diffusion coefficient. As shown in Fig. 2, this function provides a good fit of our data. In Fig. 3, we plot a fit of Eq. (8) (dashed line) to the simulation results, showing a good collapse of the bulk particle density outlines for different times. The small deviations close to the aggregate are due to the interface roughness. The ratio $`D/u`$ is called diffusion length. In our fits, this length turns out to be about $`15`$ lattice spacings, that is, about $`2`$ or $`3`$ times the lateral width of the branches for the chosen parameters. This result provides another of check the physical validity of our model, as we can compare the length obtained from the fit with that taken from Ref. . In this paper, the diffussion length is of order $`0.5`$ $`mm`$, about two times the typical branch lateral width (of order $`1`$ $`mm`$), so we may conclude that the diffusion length obtained from our model is physically consistent.

The inset in Fig. 3 shows the mean concentration front position $`z_0`$ vs. time, demonstrating that, in the stationary state, MBDLA leads to a constant velocity of the advancing front as in the experiments.

#### 2 Physical Relevance of SD

The previous subsection shows that MBDLA without SD successfully reproduces some ECD experiments, in particular, under galvanostatic conditions with not very small electric current density. However, within the MBDLA model it is impossible to understand the unexpected compactification of aggregates in low voltage experiments  or the columnar-like growth found in other situations. Unfortunately, MBDLA aggregates are always ramified at low bias. In , a phenomenological explanation of compactification was proposed by noticing that the reduction of $`s`$ leads to more compact aggregates. Therefore, it was proposed there that $`p`$ and $`s`$ should be related by a monotonous function, the simplest case being that of a linear relationship. With this procedure, reducing the bias leads to a corresponding decrease in the sticking probability, and hence to compact aggregates at low bias. However, this is an ad hoc assumption that cannot be experimentally tested, whereas its theoretical justification is not very clear. Besides that, this approximation does not reproduce other morphologies, as that reported by López-Salvans et al.  or Kahanda et al. . In view of this, it became increasingly clear that there was some crucial ingredient missing in MBDLA, and the most obvious candidate was of course SD.

At this point, it is instructive to consider carefully the work by Kahanda et al. . According to their results, as the absolute value of the overpotential decreases the aggregate becomes denser, and it is formed by several columns which are thicker at the top than at the bottom. We interpret this as a hint on the relevance of SD: If, when a particle arrives at the top of a column, it diffuses along the aggregate interface, and if the diffusion length is shorter than the column perimeter, the particle will not reach the base of the pillar or another column, with the result of a characteristic inverted triangle structure. The onset of similar triangle structures has also been reported by Pastor and Rubio . We thus came to the conclusion that it was necessary to include SD in MBDLA in order to shed further light on the nontrivial coupling of the different transport mechanisms.

#### 3 Implementation of SD in MBDLA

We have implemented SD in MBDLA in three different ways, all of them starting when a particle in the bulk (the electrolytic solution) sticks to the aggregate but has just one neighbor. We have first tried two simple irreversible rules (other similar rules yield equivalent results, so we do not include them here for brevity), named rules A and B, and an irreversible one, named rule C:

Rule A: The newly incorporated particle jumps always in the same direction, either left or right parallel to the cathode, until it reaches a site with at least $`2`$ neighbors or completes $`l`$ jumps. This rule is similar, but not identical, to the one studied in  for ballistic deposition with surface diffusion.

Rule B: In this second rule, we allow the particle to perform a random walk over the aggregate surface until it increases its coordination number, with a constant probability $`\lambda `$ to be permanently stuck to its current position (this is the so called mortal random walker ).

The last rule is characterized by Arrhenius-like jump probabilities and, what is more important, by simultaneous bulk diffusion and SD:

Rule C: This rule allows several particles to diffuse simultaneously. When a particle arrives to a coordination $`1`$ site, it sticks and jumps to one of its two nearest neighbors on the aggregate with probability $`p_n=\mathrm{exp}[E_0+(n1)E_a]`$, where $`E_{0,a}`$ are adimensional activation energies, and $`n`$ the coordination number of the target position. If the particle new position has $`2`$ or $`3`$ neighbors, it attaches to the aggregate irreversibly. Otherwise, we label the particle as a SD particle, and we allow it to take further steps. Thus, we have two kinds of diffusing particles: Particles in the bulk, distributed homogeneously with concentration $`c`$; and particles that diffuse over the aggregate surface. With probability $`r`$ we choose a bulk particle which evolves with its characteristic rules, and with probability $`1r`$ a particle on the surface which jumps to one of its nearest neighbors as we have just described for the first jump. This rule is close in spirit to the collective diffusion rules employed in studies of kinetic roughening in molecular beam epitaxy (MBE) , and in particular to MBE models beyond the solid-on-solid approximation .

The main difference between rules A and B with respect to rule C is that the latter introduces a characteristic time scale $`\tau _d=r^1`$, while in the other cases diffusion is instantaneous; then, the diffusing particle is not affected by the overall particle dynamics. As we will show below, Rule C is the only one which actually reproduces the influence of SD on the aggregates scaling and morphology. In this respect, it is important to advance that we have found that Arrhenius-like probabilities by itself are not enough to model SD: Variants of rule C with those probabilities and without the characteristic time, i.e., SD kept instantaneous, lead once again to results similar to those of rule B. All the results presented were obtained with $`E_0=3`$ and $`E_a=1`$. We have chosen these values to have jump probabilities smaller than $`1`$, but other sets of parameters yield similar results which we omit for brevity. Finally, another interesting point is that the probabilities in Rule C allow to trivially introduce temperature in the model by simply identifying $`E_{0,a}E_{0,a}^{}/k_BT`$.

Rule A, by definition, introduces a diffusion length $`l`$, but if $`l1`$ the particle jumps practically always lead to an increment of its coordination, as may be seen in Fig. 4, where some morphologies are shown for different values of $`l`$. The inverted triangle structure typical of the experiments by Kahanda et al.  is reproduced with this simple rule. Nevertheless, the tops of the pillars are unrealistically flat; another problem is that decreasing $`p`$ does not lead yet to a compact aggregate regime. Rule A is therefore not appropriate. In the case of Rule B, the diffusion length is introduced indirectly by means of the attachment probability $`\lambda `$ (see Appendix for details). The mean diffusion length can be shown to be given by $`l_D=1/(2\lambda ^{1/2})`$. Morphologies obtained with this rule are plotted in Fig. 5. Once again, and in spite of the fact that rule B allows the particles to diffuse randomly over the aggregate, the columns developed during the growth turn out completely flat at the top, and the option of rule B was excluded as well.

These pitfalls (and similar ones found by using Arrhenius-like jump probabilities, which we skip for brevity) led us to the conclusion that instantaneous SD (or limited mobility rules, in the terminology of ) is a too drastic one approximation for ECD. Taking into account that the main distinctive feature of MBDLA is its non-local character, interactions of particles diffusing along the surface with newly deposited particles are expected to be relevant. Guided by these ideas, we propose rule C, which incorporates this coupling by introducing time scales for both bulk diffusion and SD. A sample of the aggregates generated by MBDLA with rule C is shown in Fig. 6. The difference with the other two rules is immediately apparent from the plot: This more realistic rule does induce the creation of pillars as we pointed out above, this time similar to those reported by Kahanda et al.  and Pastor and Rubio  which are rough at the top. Moreover, the compactification of the aggregates at low currents appears naturally, as can be noticed by following the sequence of aggregates appearing on the same row (same value of $`r`$): Decreasing the current leads initially to less dense aggregates, until further reduction of the current gives rise to more compact aggregates. Remarkably, there is no need to change the sticking probability by hand as in MBDLA without SD or with Rules A and B. This allows us to eliminate one model parameter, the sticking probability, which we take to be $`s=1`$ from now on.

So far, we have seen that, while simple SD rules provide good results in some solid-on-solid simulation models, the complex dynamics of Laplacian systems does not allow the particles to instantaneously diffuse; rather, we must allow several particles to interact before they become permanently stuck to the aggregate. Roughly speaking, the flux of particles arriving at the aggregate defines a characteristic time $`\tau _p`$ (typically inversely proportional to the flux, i.e., to $`p`$). Once the particles have arrived to the aggregate, they diffuse until they reach a site with coordination larger than one, or equivalently, until the particle meets another diffusing particle, thus forming a dimer on the interface which cannot move anymore. A large flux of particles arriving at the interface (large $`p`$) will increase the probability of formation of those dimers, and the particles can hardly diffuse. The situation is not so simple when $`p`$ is small. On one hand, the deposition mean time $`\tau _p`$ is large, but on the other hand, the particles hardly experience the applied electric current, so the probability of attachment to a column wall before getting to the bottom of the aggregate increases. Thus the Laplacian instability is amplified leading to a compact structure formed by columns and grooves. This kind of instability has been observed in low current galvanostatic experiments . The aggregate is therefore denser but if the diffusion time is not long enough the interface is unstable. It is remarkable that this simple picture in terms of time scales allows to understand the relevance of SD in ECD experiments.

A final important remark we would like to mention is that, when the diffusion probability $`r`$ is about $`0.99`$, we have observed some evidences of what could be a morphological transition (and the subsequent change in the branches) similar to those reported by López-Salvans et al. . However, as we want to concentrate in this paper on MBDLA with SD as a generic model for all regimes of ECD experiments, we postpone a more careful study of this possibility to future work, where we will pursue the appearance of this phenomenon for different model parameters (such as $`p`$ or $`r`$).

#### 4 Electrolyte concentration

To conclude the analysis of MBDLA parameters, we show the effect of the electrolyte concentration, $`c`$. Figure 7 exhibits the morphological changes in patterns with different $`c`$ values ranging from $`0.01`$ to $`0.1`$ for different bias $`p`$ without SD. Note that when $`c0`$ the low current limit is exactly the DLA growth model . Therefore, we should keep a finite value of $`c`$ in order to diminish the unavoidable DLA characteristic instability. The results contained in the figure allow us to conclude that, insofar $`c`$ is not very small, the morphologies obtained with MBDLA do not depend strongly on the concentration, and therefore the fact that there is no direct correspondence between physical and simulated concentrations is not a drawback of the model.

### B Dynamic Scaling

The previous subsection shows that the inspection of the morphologies is a valuable method to check the validity and relevance of the model rules. Indeed, the unrealistically flat aggregates obtained with diffusion rules A and B unable them and motivate the investigation of the more realistic, non-instantaneous rule C for SD. However, in order to exploit the main virtues of MBDLA with SD and to compare with other relevant models and experiments, we must take some quantitative criteria, for example, the analysis of the interface surface roughening. To this end, let us define some functions related to the height of the aggregate at spatial position $`x`$ at time $`t`$, given by the scalar field $`h(x,t)`$. We will also review their basic features before discussing MBDLA properties.

The global width (or roughness) $`W(L,t)`$ is nothing but the rms fluctuations of the height variable $`h(x,t)`$ around its mean value $`\overline{h}_L(t)=(1/L)_xh(x,t)`$:

$$W^2(L,t)=\frac{1}{L}\underset{x}{}\left[h(x,t)\overline{h}_L(t)\right]^2,$$
(9)
where angular brackets stand for ensemble average.

Generally speaking, in many growth models, starting from $`h(x,0)=0`$ the width satisfies the dynamic scaling hypothesis of Family-Vicsek :

$$W(l,t)\{\begin{array}{ccc}t^\beta \hfill & \mathrm{if}\hfill & tL^z,\hfill \\ L^\alpha \hfill & \mathrm{if}\hfill & tL^z.\hfill \end{array}$$
(10)
The roughness exponent $`\alpha `$, the dynamic exponent $`z`$, and their ratio (growth exponent) $`\beta =\alpha /z`$, identify the universality class the model belongs to.

In the study of kinetic roughening the height-height correlation function is frequently used :

$$C^2(l,t)=\frac{1}{L}\underset{x}{}\left(h(x+l,t)h(x,t)\right)^2,$$
(11)
where ,

$$C(l,t)\{\begin{array}{ccc}t^\beta \hfill & \mathrm{if}\hfill & tl^z,\hfill \\ t^{(\alpha \alpha _{loc})/z}l^{\alpha _{loc}}\hfill & \mathrm{if}\hfill & tl^z,\hfill \end{array}$$
(12)
where $`\alpha _{loc}`$ is the so called local roughness exponent. Another important function related to the height variable $`h`$ is the power spectrum:

$$S(k,t)=\widehat{h}(k,t)\widehat{h}(k,t),$$
(13)
where $`\widehat{h}(k,t)=L^{1/2}_x[h(x,t)\overline{h}_L(t)]\mathrm{exp}(\mathrm{i}kx)`$. $`S(k,t)`$ displays a behavior consistent with the scaling form

$$S(k,t)=k^{(2\alpha +1)}s(kt^{1/z}),$$
(15)
where

$$s(u)=\{\begin{array}{ccc}u^{2\theta }\hfill & \mathrm{if}\hfill & u1,\hfill \\ u^{2\alpha +1}\hfill & \mathrm{if}\hfill & u1.\hfill \end{array}$$
(16)
The exponent $`\theta `$ takes different values depending on the type of scaling exhibited by the model. For instance, for the so called intrinsic anomalous scaling  we have $`\theta =\alpha \alpha _{loc}`$, whereas $`\theta 0`$ for Family-Vicsek scaling (including super-roughening, i.e., $`\alpha 1`$). Note that this implies $`\alpha =\alpha _{loc}`$.

To apply these ideas to MBDLA characterization, a few remarks are in order. Although, in some cases, MBDLA develops ramified aggregates leading to multivalued interfaces, i.e., interfaces with overhangs, it has been demonstrated  that the interface of the active zone in DLA simulations (the aggregate sites with larger probability of arrival) corresponds to that constructed by taking the topmost site $`h(x,t)`$ at every horizontal position $`x`$. This construction does not ensure that the measured exponents are free of interpretations , but the exponents are consistent with theoretical and experimental data . The reduction of the sticking probability $`s`$ yields denser aggregates, and overhangs do not appear at any stage of the simulation for low $`s`$ values. Besides that, if SD is present the aggregates are also more compact. In all these cases the function $`h(x,t)`$ is identical to the aggregate outline and consequently the results do not have any interpretation problem.

The main scaling features of MBDLA without SD were already reported in . Therefore, here we will briefly summarize them to facilitate comparison with results including SD, and refer the reader to  for the details. Without SD, MBDLA displays three temporal regimes: At early times the global width $`W(L,t)`$ features $`\beta =0.5`$, this value being simply due to shot noise. This stage corresponds to times at which the lateral correlation length is of the order of the lattice spacing. After this noisy transient, short and large length scales are governed by different dynamics because the bulk Laplacian field produces nonlocal effects (screening or shadowing among branches). Consequently, the local and the global roughness exponents, $`\alpha _{loc}`$ and $`\alpha `$, are different and the interface is not self-affine. The growth exponent, $`\beta `$, is larger than that of noise ($`\beta >1/2`$) because some isolated branches begin to grow independently from each other, which can be understood as a signature of the Laplacian instability. As a consequence, the interface width grows rapidly as compared with the noise fluctuations. At later times, branches spread by lateral growth and impinge upon each other. Eventually, the system reaches an asymptotic regime characterized by the Kardar-Parisi-Zhang (KPZ) universality class  exponents ($`\alpha =1/2`$, $`\beta =1/3`$, $`z=3/2`$). The KPZ equation is the paradigmatic growth model without SD, and it is given by the stochastic partial differential equation

$$\frac{h}{t}=\nu ^2h+\frac{\lambda _0}{2}\left(h\right)^2+\eta (x,t),$$
(17)
where $`\nu `$ and $`\lambda _0`$ are constants and $`\eta (x,t)`$ is a Gaussian white noise with:

$$\eta (x,t)=0,$$
(19)
$$\eta (x,t)\eta (x^{},t^{})=2D\delta (xx^{})\delta (tt^{}).$$
(20)
As mentioned above, the definition of the interface function $`h(x,t)`$ neglecting overhangs might cast some doubts  on the validity of the exponents reported in . To confirm our results, we have measured the excess velocity produced by tilting the initial substrate and imposing helicoidal boundary conditions . The inset in Fig. 8 shows that this mean velocity is well fitted by a parabola, as expected for KPZ behavior. It is important to stress that identical results are obtained using the jump rules in .

Interestingly, Schilardi et al.  report experiments with large currents (equivalent to the large values of the bias $`p`$) in excellent agreement with our model. They observe the same three time regimes: An initial transient with a behavior which could not be measured due to the resolution of the experimental device; a second transient with $`\beta >1`$ characterized by the growth of isolated branches, and a third asymptotic regime at which the interface is characterized by KPZ exponents. A plot of the mean interface velocity vs. time is also given, showing a crossover from the unstable regime to the stable one in accordance with MBDLA predictions as can be seen in Fig. 8. The global width crosses over from the instability ($`\beta >0.5`$) to $`W(L,t)t^{1/3}`$ at the time pointed out by the arrow. Note that MBDLA cannot yield $`\beta `$ larger than one, because of its discrete growth rules. This would mean that the interface width grows faster than the interface mean height. Finally, the evolution of the morphology during the experiment is also the same in MBDLA and in the experiment, as seen by comparing Fig. 9, taken from , and Fig. 10, obtained in our simulations.

We now consider the scaling behavior of MBDLA+SD for the different diffusion rules. As we have pointed out in the preceding section, large values of the diffusion length $`l`$ (rule A) generate flat aggregates. This means that $`\beta 0`$ as $`l\mathrm{}`$ at early times. Fig. 11 shows the lack of universality in the growth exponent $`\beta `$: It can be seen in this plot that $`\beta `$ decreases with $`l`$ as we expected. The same happens with rule B: As in the case of rule A, the growth exponent $`\beta `$ depends strongly on the attachment probability, $`\lambda `$ (rule B). As depicted in Fig. 11, the dependence is similar to that of model A since the diffusion length is proportional to $`\lambda ^1`$.

The scaling behavior in MBDLA with SD given by rule C is more complicated. We can recognize three different kinds of behavior, which we summarize as follows:

$`\mathbf{0.05}𝐫\mathbf{0.25}`$. The characteristic diffusion time is long, and particles diffuse rather fast along the surface (let us recall that they are picked with probability $`1r`$ at every Monte Carlo trial) without much interaction with particles arriving from the bulk, thus yielding compact aggregates, except if $`p0.05`$, because then the Laplacian field creates pillars and grooves. After a short transient the global width grows slowly and, independently of the applied current, the roughness exponents are compatible with those of the Edwards-Wilkinson universality class, whose defining equation is :

$$\frac{h}{t}=\nu ^2h+\eta (x,t).$$
(21)
Figure 12 shows the global width collapse obtained by rescaling the simulation time. The plot not only shows the Edwards-Wilkinson growth exponent, but also the $`r`$ independence of the results on a wide range of simulation parameters. Note that the collapsing time step is the one defined for MBDLA without SD divided by the characteristic diffusion time $`\tau _d=r^1`$. Figure 13 shows the collapsed power spectrum using $`\alpha =1/2`$ and $`z=2`$ (and consequently $`\beta =1/4`$) consistent with (13) with $`\theta =0`$ for Edwards-Wilkinson exponents. It is important to note that this kind of dynamic scaling has been observed in two-dimensional ECD experiments . Finally, we have to mention that the restriction $`r>0.05`$ is only due to the extremely long computational times needed to study the model for such small values of $`r`$.

$`\mathbf{0.3}𝐫\mathbf{0.7}`$. For large $`p`$, the interface is compact and grows with constant velocity. The scaling is similar to that of the preceding case. When $`p0`$, initially the interface is rough and the growth exponent $`\beta `$ is in the range $`0.350.40`$ (see Fig. 14). Some experiments have reported similar interfaces at early stages of growth : Specifically, they obtained exponents consistent with the linear MBE growth model universality class ($`\alpha =3/2`$, $`\beta =3/8=0.375`$ and $`z=4`$) that is, their interfaces could be described by the equation :

$$\frac{h}{t}=K^4h+\eta (x,t).$$
(22)
Note that for this model $`\alpha >1`$, so the interfaces generated with Eq. (22) are super-rough. In our case this short regime ceases when the mean interface height, $`\overline{h}(t)`$, is about $`8`$ to $`10`$ monolayers and the global width, $`W(L,t)`$, is about $`1`$. This is compatible with the referred experiments except that we do not observe the super-rough power spectrum. Actually, in our case the tail of $`S(k)`$ presents a time shift at large wave vectors (Fig. 15) which is incompatible  with the behavior obtained for Eq. (22). However, the basic phenomena, such as the value of the effective $`\beta `$ and the onset of the instability, are in good agreement with the experiments. After this transient, the aggregates are still compact and develop some grooves (see Fig. 6). When these grooves appear, the growth exponent $`\beta `$ rises dramatically due to the large slopes produced between grooves. Figure 16 summarizes all this by showing the variation of $`\beta `$ with time.

$`\mathbf{0.85}𝐫`$. Finally, when the diffusion time is short, three completely different situations are found as a function of the current $`p`$. For very large $`p`$, cations become ballistically driven to the aggregate and the unstable transient tends to dissappear (in fact, the $`p\mathrm{}`$ limit is the Ballistic Depostion discrete model, which is well known to belong to the KPZ universality class ). When $`p1`$ the aggregate grows as MBDLA without SD with similar parameters, except that in this case the aggregate mean density rises. That is, we succesively detect a noisy initial transient, the instability associated with the growing branches, and the KPZ asymptotic limit due to the lateral growth of the branches. The interfaces within the unstable regime (an example of which is shown in Fig. 17) are not self-affine but present intrinsic anomalous scaling : Figure 18 shows the power spectrum for $`r=0.85`$, $`p=4`$ and $`s=1`$. Figure 19 shows the collapse of the power spectrum and Fig. 20 the collapse of the height-height correlation function, $`C(l,t)`$, achieved in both cases for $`\alpha =1.78`$, $`\alpha _{loc}=0.49`$, $`z=2.51`$ and $`\beta =0.71`$.

For intermediate $`p`$ values (between $`0.25`$ and $`1`$, for almost every $`r`$), the aggregate is formed by several compact thin branches which grow vertically and parallel to one another. In this case, the notions of a rough surface or dynamic scaling are meaningless. Finally, for small $`p`$ some compact branches grow at the expenses of the others, so typically, one or two branches grow more than the others. As in the case of the parallel branches, it is meaningless to talk about interface roughening.

## IV Discussion and Conclusions

Our first conclusion is that MBDLA is a simple computational model which incorporates in a natural way some of the basic mechanisms involved in ECD experiments. The original model  was already known to be in good agreement with some experiments . In this paper, we have provided much more evidence showing that MBDLA explains some of the morphological changes due mainly to the applied electric current and, what is more important, it predicts the recently observed KPZ scaling behavior in the high current limit (for which SD is not too relevant)  and observed also at low currents . Before this regime is reached, there is an unstable transient within which MBDLA interfaces present intrinsic anomalous scaling. We believe this type of scaling is due to SD not being able to communicate different portions of the interface fast enough, so that they grow independently from one another. This is analogous of the anomalous scaling occurring in the non-linear surface diffusion equation studied in . In our case, the different portions feature a value of the roughness exponent $`\alpha _{loc}0.5`$, similarly to the interface subject to columnar disorder studied in .

Secondly, the main point of our paper is that, as we have seen, MBDLA without SD cannot explain low current experiments in which the characteristic dense branching aggregates of high current experiments are replaced by compact and column-like aggregates. Our working hypothesis was that the latter kind of patterns are due to the competition between the Laplacian field of the cations in the disolution and the SD current on the interface. Thus, our ECD model, which we wanted to improve as to explain, at least qualitatively, the complete ECD phenomenology, should incorporate a new rule for the diffusion of the adatoms attached to the aggregate. Hence, we have tried out some SD rules similar to those often used in growth models for molecular beam epitaxy . We have verified that instantaneous diffusion rules, namely rules that “freeze” the bulk particles while the most recently attached particle finds its way through the surface, do not lead to correct results in the low current limit, and produce very unrealistic, flat-topped morphologies. We have then been forced to conclude that the nonlocal character of MBDLA demands a diffusion rule which couples the overall cation dynamics: This is the rule we have named C. It introduces a characteristic diffusion time $`\tau _d=r^1`$ which competes with the time scale related to the net flux of particles arriving to the interface (which, in fact, is proportional to the applied electric current density). With this SD rule, the morphologies at low, medium and high currents are compatible with those observed by Trigueros et al.  for low, medium and high applied voltages respectively. This diffusion time $`\tau _d`$ cannot be controlled from the experimental point of view, but fortunately there are wide ranges of parameters over which the simulated morphologies hardly change, which means that the description of the experiments provided by MBDLA with SD is robust and does not need an incontrollable parameter to be tuned. We have also compared MBDLA with SD with the experiments reported by Pastor and Rubio , which characterize the product interfaces by the MBE exponents. MBDLA seems to reproduce the latter behavior for very short times and short length scales as can be seen in Figs. 14 and 15, but these results are not too significative, as they are not as accurate as we would need to make any strong claim, and could be due to the appearance of a characteristic short length scale. There is a difficulty in this respect, as MBDLA scaling is intrinsically anomalous whereas the results in  support standard, super-rough scaling. With the data presently at hand we have to conclude that MBDLA with SD does not describe quantitatively all the aspects of the very low current regime, but the fact it does describe most of them and, above all, the compactification of the aggregates, makes us confident that MBDLA with SD is a very good general model for ECD.

To conclude, we note that the model presented has the basic ingredients of ECD phenomena: Diffusion, electromigration and surface diffusion, but for this reason, we have to pay a big price in terms of computational time. MBDLA without SD is a very time-consuming model, and the diffusion rules make the analysis and the simulations a patience exercise. It has certainly been an improvement to find that SD rule C allows us to skip the sticking probability parameter, thus resducing the parameter space, but even then if averages of relevant quantities over large ensembles are required, a great deal of computational resources would be needed. Of course, this disadvantage can be removed by a careful reprogramming of the algorithm, but that is another line of research. As our goal was to identify the most important factors involved in ECD, we do believe that despite the computing limitations of the model, MBDLA with SD is a powerful tool to repoduce some unclear features of this kind of growth experiments, and has helped understand what are the most relevant transport properties and how they couple in different parameter regions. We hope that this work suggests further experiments to find out whether MBDLA with SD is the complete, general model for ECD or else if there are still regions which need separate modelling.

###### Acknowledgements.

It is a pleasure to thank Miguel Angel Rodríguez and Francesc Sagués for their help through many discussions on the modeling of ECD and on the subtleties of dynamic scaling during the last few years, as well as for a critical reading of the manuscript. In addition, we are thankful to Javier Buceta, Juanma Pastor, Javier de la Rubia, and Miguel Angel Rubio for sharing with us their results and expertise on ECD at very low currents. We would also like to thank Roberto Salvarezza and Leticia Schilardi for their kind permission to reproduce Fig. 9 from . We are indebted to José Cuesta and Ricardo Brito for helpful suggestions about the simulation rules and to the members of GISC for their interest in this work. Work at GISC has been supported by DGES (Spain) Grant No. PB96-0119 and CAM Grant No. 07N/0034/98.

## Surface Diffusion Rule B

In Rule B, surface diffusion starts when a particle has first arrived to the aggregate and has attached to it (with probability $`s`$). The particle jumps with equal probability to one of its two nearest neighboring sites on the aggregate until it increases its coordination number. The particle has an additional probability, $`\lambda `$, of being permanently attached. This kind of particle is usually termed a mortal random walker . The random walk is performed between two absorbing boundaries, namely, a couple of sites with higher coordination (two or three). One could try to determine a priori the total number $`N`$ of jumps the particle has to perform in each realization drawing such number from the probability for the particle to take $`N`$ steps on a flat line if it avoids to stick $`N1`$ times and die at the $`N`$th jump. This probability can be easily calculated to be given by

$$P_N(\lambda )=\lambda (1\lambda )^{N1}.$$
(23)
However, the absorbing boundaries disallow this procedure. In any case, we have compared the simulation results by allowing the particle to perform an actual mortal random walk, and to perform a simple random of $`N`$ steps given by Eq. (23). Both results are hardly different. Thus, we can approximately calculate from Eq. (23) the mean and variance of the maximum number of jumps, given by:

$`\overline{N}={\displaystyle \frac{1}{\lambda }},`$ (24)
$`\sigma _N={\displaystyle \frac{\sqrt{1\lambda }}{\lambda }}.`$ (25)
For a flat interface, the particle mean position would be $`0`$ but its variance would be

$$\sigma _N=\overline{N}^{1/2}/2=1/(2\lambda ^{1/2}),$$
(26)
which provides the characteristic diffusion length $`l_D=1/(2\lambda ^{1/2})`$.


format: text

text:
# References

MAGNETIC BEHAVIOR OF A MIXED ISING FERRIMAGNETIC MODEL IN AN OSCILLATING MAGNETIC FIELD.
G. M. Buendía, E. Machado

Departamento de Física. Universidad Simón Bolívar,
Apartado 89000, Caracas 1080, Venezuela.

ABSTRACT

The magnetic behavior of a mixed Ising ferrimagnetic system on a square lattice, in which the two interpenetrating square sublattices have spins $`\sigma `$ ($`\pm 1/2`$) and spins $`S`$ ($`\pm 1,0`$), in the presence of an oscillating magnetic field has been studied with Monte Carlo techniques. The model includes nearest and next-nearest neighbor interactions, a crystal field and the oscillating external field. By studying the hysteretic response of this model to an oscillating field we found that it qualitatively reproduces the increasing of the coercive field at the compensation temperature observed in real ferrimagnets, a crucial feature for magneto-optical applications. This behavior is basically independent of the frequency of the field and the size of the system. The magnetic response of the system is related to a dynamical transition from a paramagnetic to a ferromagnetic phase and to the different temperature dependence of the relaxation times of both sublattices.

INTRODUCTION

The behavior of ferrimagnetic compounds in the presence of oscillatory fields have long been used for technological applications such as high-density magneto-optical recording , but little is known about the mechanisms responsible for this behavior. In a ferrimagnet the different temperature dependencies of the sublattice magnetizations raise the possibility of the appearance of compensation temperatures: temperatures below the critical point, where the total magnetization is zero . It has been shown experimentally that the coercive field is very strong at the compensation point favoring the creation of small, stable, magnetic domains . This temperature dependence of the coercivity near the compensation point can be applied to writing and erasing in high-density magneto-optical recording media, where the temperature changes are achieved by local heating the films by a focused laser beam. It has been shown that magneto-optic thin films with compensation temperatures higher than room temperatures can attain a direct overwrite capability . As far as we know there has been only very few crude attempts to reproduce theoretically the increase of the coercivity near the compensation point using mean-field approaches .Recently new classes of magnets are being synthesized with molecular organic chemistry techniques . Biocompatible, organic materials, optically transparent, with spontaneous moments at room temperature are not far from reality. Ferrimagnetic ordering seems to play a fundamental role in some of these materials. Ferrimagnetic compounds called Prussian blue analogs, with a critical temperature of 240 K have been reported . Organometallic compounds as the amorphous $`V(TCNE)_xy(solvent)`$ where TCNE is tetracyanoethylene are believed to have ferrimagnetic structure and ordering temperatures as high as 400 . Some of these compounds have compensation temperatures near 30 K . Most of these compounds have been synthesized by assembling molecular building blocks of different magnetic moments in such a way that adjacent magnetic moment are antiparallel . Since real ferrimagnets have extremelly complicated structures mixed Ising models have been introduced as simple systems that can show ferrimagnetic behavior  and may show compensation points when their Hamiltonian includes second-neighbor interactions . In this article, we present a Monte Carlo study of a mixed Ising spin system, where spins that can take the values $`\pm 1/2`$, and spins that can take the values $`\pm 1,0`$, are nearest neighbors on a 2-dimensional square lattice and interact antiferromagnetically. Spins of the same type are next-nearest neighbors. We analyze the magnetic response of this system in the presence of an oscillating magnetic field. From these studies, we determine the dynamic order parameter, the coercive field and their variation with temperature, frequency and amplitude of the applied field, and size of the system. The results reproduce the rapid increase of the coercitivity at the compensation temperature. The dynamical order parameter calculations suggest that the model exhibits a phase transition between a paramagnetic and a ferromagnetic region. A similar result was observed by a mean-field study of a simpler version of this model . Mean-field approaches and Monte Carlo simulations indicates the presence of a dynamical phase transition in a kinetic Ising model . However the distinctive behavior of the coercive field at the compensation temperature in ferrimagnets seems to be related to the different relaxation times of the sublattices.

THE MIXED ISING MODEL

Our model consists of two interpenetrating square sublattices. One sublattice has spins $`\sigma `$ that can take two values $`\pm 1/2`$, the other sublattice has spins $`S`$ that can take three values, $`\pm 1,0`$. Each $`S`$ spin has only $`\sigma `$ spins as nearest neighbors and vice versa.

The Hamiltonian of the model is given by,

$$=J_1\underset{nn}{}\sigma _iS_jJ_2\underset{nnn}{}\sigma _i\sigma _k+D\underset{j}{}S_j^2H(t)\left(\underset{i}{}\sigma _i+\underset{j}{}S_j\right)$$
(1)
where the $`J`$’s are exchange interaction parameters, $`D`$ is the crystal field and $`H`$ is an oscillating magnetic field of the form,

$$H(t)=H_0cos(\omega t)$$
(2)
where $`\omega `$ is the frequency of the external field, its period is given by, $`\mathrm{\Theta }=2\pi /\omega `$. The $`J`$’s, $`D`$ and $`H_0`$ are all in energy units. We choose $`J_1`$$`=`$$`1`$ such that the coupling between nearest neighbors is antiferromagnetic.

Previous results with Monte Carlo and Transfer-Matrix techniques have shown that the $`J_1`$$``$$`D`$ model ($`J_2`$ and $`H`$ are equal to zero) does not have a compensation temperature. These studies show that a compensation temperature is induced by the presence of the next-nearest-neighbor (nnn) ferromagnetic interaction, $`J_2`$, between the $`\pm 1/2`$ spins. The minimum strength of the $`J_2>0`$ interaction for a compensation point to appear depends on the other parameters of the Hamiltonian .

MONTE CARLO CALCULATIONS

We use standard importance sampling techniques to simulate the model described by Eq. (1) on a $`L\times L`$ square lattice with periodic boundary conditions. Configurations are generated by randomly choosing spins on the lattice and flipping them one at a time according to a heat bath algorithm. In each complete sweep through the lattice $`L\times L`$ sites are visited. Each Monte Carlo step per spin is associated to a time interval, $`\tau _S`$, such that the frequency of the external field can be written as,

$$\omega =\frac{2\pi }{(NMCS)\tau _S}$$
(3)
where $`NMCS`$ is the number of Monte Carlo steps per spin necessary to cover an entire cycle of the field. To perform the simulations we arbitrarily choose $`\tau _S`$ to be one, such that $`\mathrm{\Theta }=NMCS`$. Our program calculates the sublattice magnetizations per site at the time $`t`$ defined as,

$$M_1(t)=\frac{2}{L^2}\underset{j}{}S_j(t),M_2(t)=\frac{2}{L^2}\underset{i}{}\sigma _i(t)$$
(4)
and the total magnetization per spin at the time $`t`$, $`M(t)=\frac{1}{2}\left[M_1(t)+M_2(t)\right]`$. The averages are taken over all configurations, the sums over $`j`$ are over all sites with $`S`$ spins, and the sums over $`i`$ are over all sites with $`\sigma `$ spins. Each sum has $`L^2/2`$ terms.

The compensation temperature is defined as the temperature below the critical, $`T_{\mathrm{comp}}<T_{\mathrm{crit}}`$, where the two sublattice magnetizations cancel each other such that the total magnetization is zero, i.e.,

$$|M_1(T_{\mathrm{comp}})|=|M_2(T_{\mathrm{comp}})|$$
(5)
and

$$\mathrm{sign}[M_1(T_{\mathrm{comp}})]=\mathrm{sign}[M_2(T_{\mathrm{comp}})].$$
(6)
To characterize the time behavior we calculate the dynamical order parameter $`Q`$ defined as,

$$Q=\frac{2\pi }{\omega }M(t)𝑑t.$$
(7)
The closed integral implies that the integral is performed over a cycle of the external magnetic field.

RESULTS

The value of Q is calculated by averaging its values over 100 cycles of the external field, once the system is in its stationary state. Most of the measurements were done for a $`L=40`$ lattice. Lattices of different sizes were used to study the finite-size effects. In Fig. 1 we show a hysteresis loop, $`M(t)`$ vs $`H(t)`$, for a particular combination of parameters in the Hamiltonian. The coercive field $`H_c`$ is defined as the minimum value of the external field needed for the total magnetization to go to zero, as is indicated in the figure. In Fig. 2 and Fig. 3 we show the coercive field vs the temperature for oscillating fields of several amplitudes, $`H_0`$. In the same figures we also plot the total magnetization for the equivalent system subject to a constant field of magnitude $`H_0`$. Notice that the compensation temperature, defined as the point where the total magnetization is zero, previous verification that Eq. 5 and 6 are satisfied, increases with the magnetic field, whereas the temperature at which the magnetization becomes discontinuous does the opposite. At a certain field, which amplitude depends on the parameters of the Hamiltonian, both temperatures become equal and for any field of larger amplitude there is no more compensation point, as can be seen in Fig. 3(b). From the figures it is clear that the coercive field increases in the vicinity of the compensation temperature where it reaches its maximum . These results are summarized in Fig. 4. As expected, the maximum value of the coercive field at the compensation temperature is given by $`H_0`$.

It is interesting to notice the asymmetric behavior of the coercive field around the compensation point. In the low temperature region, $`T<T_{\mathrm{comp}}`$, the coercive field decreases with increasing $`T`$ until it reaches a minimum, after which grows rapidly reaching its maximum at $`T_{\mathrm{comp}}`$, when $`T>T_{\mathrm{comp}}`$ the coercive field decreases. Notice that for small values of $`H_0`$ there is a range of temperatures for which the coercive field is not defined (see Fig. 4). This behavior of the coercive field has been observed experimentally . This result can be understood by looking at Fig. 5 were it is shown how the hysteresis loop changes with the temperature. As the temperature increases the loop moves in such a way that the coercive field increases until it reaches its maximum, after which, if the temperature keeps increasing, the loop stays below (or above) the $`M=0`$ axis without crossing it, meaning that the applied field is not strong enough to flip the spins. If we look at Fig.6 where we plot the coercive field and the dynamical order parameter vs the temperature we see that there is a dynamical phase transition between a paramagnetic region, $`Q0`$ and a ferromagnetic region $`Q0`$, the region where the coercive field is not defined is well into the ferromagnetic phase where the magnetization does not changes sign.

By changing the size of the system and studying the behavior of the coercive field, see Fig. 7, we notice that there are finite size effects, particularly evident for small systems ($`L<20`$). However, for larger systems, the location of the peak of the coercive field around the compensation temperature seems to be independent on the size of the system. For small systems ($`L<20`$) the peak of the coercive field appears before the system reaches its compensation temperature. Also the coercive field seems to decrease more rapidly for the larger systems.

In Fig. 8 we present some results that show the dependence of the coercive field with the size of the system, these results agree qualitatively with the experimental behavior of magnetic films and nanostructured Fe and Ni samples  for which the coercitivity depends on the average size of the grain. The size dependence of the coercive field is very similar to the size dependence of the switching field of a kinetic Ising model (field at which magnetization reversal is thermally induced on experimental timescales for given temperatures and system sizes), which behavior has been shown to be strongly dependent on the modes by which the system decays .

Next, we explore how the results depend on the frequency. In Fig. 9 we show how the coercive field vs the temperature changes for different values of the frequency of the external field. We found a quite different response to the frequency of the magnetic field depending on the dynamical phase of the system. In the paramagnetic phase ($`Q0`$, see Fig. 6) the coercive field is larger for systems driven by fields with higher frequency, but in the ferromagnetic phase ($`Q0`$) just the opposite happens, as can been seen in Fig. 9. This behavior is related to the temperature dependence of the relaxation time in the different regions . In the ferromagnetic phase we must take into account the relaxation time of both sublattices the $`\sigma `$ and the $`S`$, whereas in the paramagnetic phase only the $`\sigma `$ one is relevant because the $`S`$ lattice follows the field with almost no delay . We also notice in Fig.9 that when the field has a high frequency the maximum value of the coercive field (that occurs at the compensation temperature) does not reach the field amplitude, i.e., the coercive field does not reach its saturation value.

In Fig. 10 we show the behavior of the coercive field vs the inverse frequency for different values of the temperature. If the field has a long period the coercive field seems to reach a value that is independent of the frequency and depends on the temperature. Again we expect that this behavior is related to the temperature dependence of the relaxation times in the different phases.

CONCLUSIONS

We have applied a Monte Carlo algorithm to the study of the magnetic response of a mixed Ising ferrimagnetic model to an oscillating magnetic field. We found that this model gives very good qualitatively agreement with the magnetic behavior of real ferrimagnets. It shows a rapid increase of the coercive field at the compensation temperature, a crucial feature that makes ferrimagnetic compounds extremely useful for thermo-optical applications. It also reproduces qualitatively the dependence of the coercitivity with the size of the sample observed experimentally. The results show the existence of a dynamical phase transition in which the mean-period averaged magnetization, $`Q`$, changes from $`Q0`$ to a $`Q0`$. Work in progress indicates that, as recent studies shows is also the case for the kinetic Ising model , some aspects of the hysteretic response as its dependence on the frequency and amplitude of the oscillating field, depends on the metastable decay mode. Explaining the different behavior in the different regimes depending on the nucleation mechanism (i.e., single-droplet or multidroplet ) by which the system decays c̃iteMachado2.

ACKNOWLEDGMENTS

We are indebted to Mark Novotny and Per An Rikvold for many insightful comments during the course of this work. G. M. B also acknowledges the kind hospitality of the Supercomputer Computations Research Institute of Florida State University at Tallahassee, Florida.

FIGURE CAPTIONS

Figure 1. Hysteresis loop ($`J_2=6`$, $`D=1.9`$, $`k_BT=0.5`$, $`\omega =\pi /30`$). The coercive field is indicated.

Figure 2. Coercive field and magnetization vs temperature. ($`J_2=6`$, $`D=1.9`$, $`\omega =\pi /100`$). a)$`H_0=0.1`$ b)$`H_0=0.5`$. Notice the discontinuity in the coercive field.

Figure 3. Coercive field and magnetization vs temperature. ($`J_2=6`$, $`D=1.9`$, $`\omega =\pi /100`$). a)$`H_0=0.8`$ b)$`H_0=2.1`$. For this choice of parameters there is not compensation point for $`H_0>1`$.

Figure 4. Coercive field vs temperature ($`J_2=6`$, $`D=1.9`$, $`\omega =\pi /100`$). The maximum value of $`H_\mathrm{c}`$ is given by $`H_0`$.

Figure 5. Hysteresis loop ($`J_2=6`$, $`D=1.9`$, $`\omega =\pi /100`$, $`H_0=0.5`$). Notice that for high temperatures there is no coercive field (see Fig. 6).

Figure 6. Coercive field and dynamical order parameter vs temperature ($`\omega =\pi /30`$, $`H_0=0.5`$, $`J_2=6`$, $`D=1.9`$).

Figure 7. Coercive field vs temperature for different lattice sizes ($`\omega =\pi /100`$, $`H_0=0.5`$, $`J_2=6`$, $`D=1.9`$).

Figure 8. Coercive field vs lattice size ($`\omega =\pi /100`$, $`H_0=0.5`$, $`J_2=6`$, $`D=1.9`$). The lines are guides for the eye.

Figure 9. Coercive field vs temperature ($`H_0=0.5`$, $`J_2=6`$, $`D=1.9`$).

Figure 10. Coercive field vs inverse frequency ($`H_0=0.5`$, $`J_2=6`$, $`D=1.9`$). The lines are guides for the eye.


text:
# 1 Introduction

## 1 Introduction

The longbase laser interferometric gravitational wave detectors are under construction at present time \[1-3\]. Their sensitivity to metric perturbation will be about $`h10^{21}`$ that corresponds to the classical regime of operation. However for future installations with projected sensitivity $`10^{22}÷10^{23}`$ the quantum features of the measurement process can play a significant role. At the same time there are no limits of principle on the accuracy of measurement of external classical force. Therefore the methods and schemes which give the possibility to overcome the quantum measurement limitations (or the so called standard quantum limit, SQL) is of vital importance for future generation of gravitational wave experiments.

There are several procedures which allow to achieve the sensitivity larger than the SQL . For example in  an optimal filtration procedure for the simplest variant of the optical sensor \- a mirror attached to a mechanical resonator and illuminated with a coherent pump field - was considered. An external force acting on the mechanical oscillator displaces its equilibrium position and thus changes the phase of the reflected field. The vacuum fluctuations of the input light act on the oscillator through the radiation pressure effect and constitute the back action noise of the measuring apparatus. For such system two quadratures of reflected wave are correlated. Using correlation (phase sensitive) processing of two quadratures one can increase a signal-to-noise ratio and overcome the SQL.

However the gain in sensitivity for the schemes overcoming the SQL is usually proportional to the square root of the ratio of laser power used for pumping the interferometer and an optimal power that corresponds to the point where the sensitivity of the interferometer achieve the SQL \[4-6\]. Unfortunately the optimal power is impracticably large, about several dozens of kilowatts that restrains the experimental implementation of the technique.

The pumping with the ultrashort periodic laser pulses can be technically advantageous over a continuous wave pumping for practical realization of the schemes overcoming the SQL. Actually for a large power a problem of generating a train of short high-intensity laser pulses can be technically easier than a problem of cw light generation (when the averaged powers for two cases are equal) because in the first case the energy in laser resonator is spread over the large frequency band (and different spatial longitudinal modes) and high intensities can be produced relatively easily. At the same time the amplitude and frequency stability of the pulsed pump in the case of a mode locked laser can be at the same level as for the monochromatic pump . For example in  the stability of intermode beats for the mode locked laser output was estimated as $`510^{12}`$ in 10 s.

Another consideration is that the perspectives of squeezed states generation with high nonclassicality seem more realistic for the case of short high-intensity laser pulses allowing the use of squeezed pulsed pump in displacement transducers .

Finally an analog to digital conversion is usually used in modern experiment during the processing of the output. Therefore it seems natural to take the pulsed pump at once so that the output will comprise a set of the values for appropriate variable at definite times.

The goals of this article are to consider a displacement transducer consisting of a mirror attached to a mechanical oscillator and illuminated with a train of high-intensity laser pulses, to reveal the algorithm of optimal signal processing for such transducer and to estimate the sensitivity of the scheme to a measurement of classical external force.

The model of displacement transducer and basic equations of motion is considered in section 2. The sensitivities for traditional measurement scheme and for correlative processing of the output quadratures in the case of time independent pump are estimated in sections 3 and 4 correspondingly. The pulsed pump for the displacement transducer is considered in section 5. The conclusions are in section 6.

## 2  Model for displacement transducer and transformation of quadrature components

Let consider the most simple case of optical displacement transducer - a mirror attached to a mass of a mechanical oscillator and illuminated with a train of high-intensity laser pulses. An external force displaces an equilibrium position of mechanical oscillator changing the phase of reflected wave. The variation of the reflected field phase is measured by a readout system. This model is easy to calculate and it contains at the same time all features of displacement transducers with pulsed pump. For the incident $`E_\mathrm{i}`$ and reflected $`E_\mathrm{r}`$ waves one can use the quasimonochromatic approximation

$`E_\mathrm{i}`$ $`=`$ $`(A(tx/c)+a_1)\mathrm{cos}\omega _\mathrm{p}(tx/c)a_2\mathrm{sin}\omega _\mathrm{p}(tx/c)`$
$`E_\mathrm{r}`$ $`=`$ $`(B(t+x/c)+b_1)\mathrm{cos}\omega _\mathrm{p}(t+x/c)b_2\mathrm{sin}\omega _\mathrm{p}(t+x/c)`$ (1)
where $`A(tx/c)`$ and $`\omega _\mathrm{p}`$ are an amplitude (mean value) and a frequency of the pump wave, $`a_1`$ and $`a_2`$ are the operators of the quadrature components (fluctuations) of the pump wave (vacuum for coherent state$`),B(t+x/c)`$ is an amplitude (mean value) of the reflected wave, $`b_1`$ and $`b_2`$ are the operators of the quadrature components (fluctuations) of the reflected wave. The periodic envelope function $`A(tx/c)`$ consists of a train of equally spaced pulses and the duration of each pulse is considerably larger than the period of light wave but considerably smaller than the period of the mechanical oscillator.

To obtain the equation coupling the amplitudes of the incident and reflected waves for the moving mirror one can use a transformation of electromagnetic field for moving reference frame . For a constant velocity of the mirror $`V`$ one has

$$E_\mathrm{r}=[(1V/c)/(1+V/c)]E_\mathrm{i}\mathrm{exp}(2\mathrm{i}\omega _\mathrm{p}X/c)$$
(2)
where for simplicity the reflection coefficient of the mirror is taken to be $`r1`$ and $`X`$ is the position of the mirror. Let suppose that this expression is valid also for the slowly varying velocity $`V(t)`$ and position $`X(t)`$ of the mirror and $`V(t)c`$ (the validity of equation (2) has been proved for the mirror consisting of free electrons for the general case of relativistic velocity $`V(t)`$ in ). Then in linear approximation in $`V/c`$ one can obtain from equation (2) the following expression

$$E_\mathrm{r}=(12V(t)/c2\mathrm{i}\omega _\mathrm{p}X(t)/c)E_\mathrm{i}$$
(3)
The first term in (3) is an amplitude modulation of the reflected wave due to the mirror movement and the second is a phase modulation. For slow motion of the mirror $`V\omega _\mu X`$ ($`\omega _\mu `$ is a frequency of mechanical oscillator) and the second term in brackets is considerably smaller than the third term. Therefore for the transformation of the quadrature components of the field one can obtain

$`b_1(t)`$ $`=`$ $`a_1(t)`$
$`b_2(t)`$ $`=`$ $`a_2(t)+2\omega _\mathrm{p}A(t)X(t)/c`$ (4)
For the equation of mirror motion one has

$$\ddot{X}(t)+2\delta _\mu \dot{X}(t)+\omega _\mu ^2X(t)=M^1(F_\mathrm{s}(t)+F_\mathrm{p}(t)+F_{\mathrm{th}}(t))$$
(5)
where $`M`$ and $`\delta _\mu `$ are the mass and the damping coefficient of mechanical oscillator, $`F_\mathrm{s}(t)`$ is a signal force, $`F_\mathrm{p}(t)`$ is radiation pressure force and $`F_{\mathrm{th}}(t)`$ is a force associated with the damping of the oscillator. Let suppose for simplisity that $`\delta _\mu `$ tends to zero. Then the displacement $`X(t)`$ of the mirror will consist of two parts - a signal displacement $`X_\mathrm{s}(t)`$ and a radiation pressure displacement $`X_\mathrm{p}(t)`$. For $`F_\mathrm{p}(t)`$ one has

$$F_\mathrm{p}(t)=SA(t)a_1(t)/(4\pi )$$
(6)
where $`S`$ is a cross section of the laser beam. Therefore the equations of motion for the displacement transducer are

$`b_1(t)`$ $`=`$ $`a_1(t)`$
$`b_2(t)`$ $`=`$ $`a_2(t)+2\omega _\mathrm{p}A(t)X(t)/c`$ (7)
$`\ddot{X}(t)`$ $`+`$ $`2\delta _\mu \dot{X}(t)+\omega _\mu ^2X(t)=M^1(F_\mathrm{s}(t)+SA(t)a_1(t)/(4\pi ))`$
## 3  Sensitivity for a traditional measurement scheme

For traditional measurement scheme  the amplitude of the pump is constant. Therefore one can easily obtain the transformation relations for the quadratures $`b_1`$ and $`b_2`$ from equations (7)

$`b_1(\omega )`$ $`=`$ $`a_1(\omega )`$
$`b_2(\omega )`$ $`=`$ $`a_2(\omega )+\lambda \xi (\omega )A^2a_1(\omega )+A\xi (\omega )F_\mathrm{s}(\omega )`$ (8)
where $`\xi (\omega )=2\omega _\mathrm{p}G(\omega )/c`$, $`G(\omega )=\left[M(\omega ^22\delta _\mu \mathrm{i}\omega +\omega _\mu ^2)\right]^1`$ is mechanical oscillator transfer function and $`\lambda =S/(4\pi )`$.

Only quadrature $`b_2`$ contains the signal and it is this quadrature that is measured in traditional measurement scheme . This corresponds to the measurement of the phase of the reflected wave. The first term in the right hand side of equation (8) for $`b_2`$ can be treated as an additive noise and the second term as a back action noise. For small pump amplitudes the sensitivity is increasing with the increase of $`A`$ because the signal is proportional to $`A`$. However for large pump amplitudes the second term in r.h.s. of (8) becomes dominant and the sensitivity is decreasing with the increase of $`A`$. Therefore there is an optimal value of the pump amplitude and the sensitivity to external force at this pump amplitude is just the SQL .

## 4 Correlative processing of quadratures for time independent amplitude of the pump

Two quadratures of the reflected field according to equation (8) have the dependence on the amplitude fluctuations of the incident field $`a_1`$. Therefore one can expect that the sensitivity can be increased for the correlative processing of the output . Actually if one combine with appropriate weight coefficients the quadratures $`b_1`$ and $`b_2`$ of the output wave then in this combination the noise term depending on $`a_1`$ can be cancelled. This weighting can be done by a homodyne detector with appropriate choise of a local oscillator phase $`\varphi `$.

Let the field of the local oscillator have the form

$$E_\mathrm{L}(t)=A_\mathrm{L}\mathrm{cos}(\omega _\mathrm{p}t+\varphi )$$
(9)
Then the photodetector output is proportional to the following expression according to (1), (9)

$$I_{\mathrm{pd}}A_\mathrm{L}(b_1\mathrm{cos}\varphi +b_2\mathrm{sin}\varphi )$$
(10)
and at certain frequency $`\omega _\mathrm{f}`$ one can obtain

$$I_{\mathrm{pd}}A_\mathrm{L}[a_1(\omega _\mathrm{f})(\mathrm{cos}\varphi +\lambda \xi (\omega _\mathrm{f})A^2\mathrm{sin}\varphi )a_2(\omega _\mathrm{f})\mathrm{sin}\varphi +A\xi (\omega _\mathrm{f})F_\mathrm{s}(\omega _\mathrm{f})\mathrm{sin}\varphi ]$$
(11)
Therefore choosing the phase $`\varphi `$ according to the equation ($`\xi (\omega _\mathrm{f})`$ is real for $`\delta _\mu =0`$)

$$\mathrm{cos}\varphi +\lambda \xi (\omega _\mathrm{f})A^2\mathrm{sin}\varphi =0$$
(12)
one can make the photocurrent insensitive to the amplitude fluctuations $`a_1`$ of the input field at certain frequency $`\omega _\mathrm{f}`$ of the signal. In this case the increase of the pump amplitude $`A`$ results in the relative increase of the output signal at frequency $`\omega _\mathrm{f}`$ according to the equation (11) with respect to the noise level defined by $`a_2`$.

For compensation of the back action noise inside definite frequency band one has to use the time dependent local oscillator phase $`\varphi (t)`$ . In this case the optimal dependence of $`\varphi `$ on $`t`$ is defined by the displacement transducer transfer function $`\xi (\omega _\mathrm{f})`$ and by the spectrum of the external force $`F_\mathrm{s}(\omega )`$ .

So a signal-to-noise ratio is proportional to $`A^2`$ (there is no optimal power) and in principle there is no sensitivity limitation by the SQL. In real experiment when the pump power gets larger the output signal and noises get smaller according to equation (11) if the condition (12) is kept valid therefore when $`A`$ becomes greater than a certain value then the noises of photodetector electronics can limit the sensitivity. However this noises have technical character and will be neglected in the following.

Another sensitivity restriction can arise due to the damping in mechanical oscillator (mirror) . This problem is general for all supersensitive measurements. At the same time an intrinsic dissipation obtained in modern experiments for mechanical oscillator is far larger (by several orders of magnitude) than the value expected from the first principles  therefore it can be treated also as a technical problem now and will not be adressed below.

It is worth to mention that the increase in sensitivity over the usual measurement scheme occurs here due to utilization of the internal squeezing (self-squeezing) of the reflected beam because of the nonlinear (quadratic) interaction of the incident light and the mirror . Actually two quadratures of the reflected beam are correlated and it is this fact that allow to use the correlative processing of the output. On the other hand the correlation of the quadratures according to equations (8) means the squeezing of the beam and the larger the correlation coefficient $`\lambda \xi (\omega )A^2`$ the larger the internal squeezing .

## 5 Sensitivity for the pulsed pump

Let consider the periodic envelope $`A(t)`$ which consists of a train of equally spaced pulses with duration $`\tau `$ and period $`T`$. The spectrum of this pump has also the form of a train of pulses in frequency domain with the distance between neighbour pulses

$$\omega _\mathrm{q}=2\pi T^1$$
(13)
For the amplitude of the pump $`A(t)`$ one can use now the expansion into the Fourier series

$$A(t)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}g_n\mathrm{exp}(\mathrm{i}n\omega _\mathrm{q}t)$$
(14)
and the particular form of $`A(t)`$ is defined by the set of Fourier amplitudes $`g_n`$.

The response of the displacement transducer now have many frequency components at $`\omega =n\omega _\mathrm{q},n=0,1\mathrm{}`$ according to the equations (4) and each frequency component contains the signal part besides the radiation pressure force $`F_\mathrm{p}(t)`$ have also wide spectrum (cf. (6)). So there are two problems: how to collect the signal parts from the whole spectral band of the output and how to achieve the compensation of the radiation pressure noise in the output. It is clear that the monochromatic local oscillator is inappropriate for the homodyning because quadrature $`b_1(t)`$ of the output signal contains in this case the quadrature $`a_1(t)`$ of the input noises only from one frequency and the radiation pressure force $`F_\mathrm{p}(t)`$ in expression for $`b_2(t)`$ (cf. (7)) contains $`a_1(t)`$ from all frequencies $`n\omega _\mathrm{q}`$ therefore the full compensation is impossible.

Fortunately two problems can be overcome by the use of the pulsed local oscillator with the amplitude time dependence resembling that for the pump.

For the radiation pressure displacement $`X_\mathrm{p}`$ of the mechanical oscillator one has from equations (5), (6) and (14) the following expression

$$X_\mathrm{p}(\omega )=G(\omega )F_\mathrm{p}(\omega )=\lambda G(\omega )\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}g_na_1(\omega n\omega _\mathrm{q})$$
(15)
For the quadrature transformation one can obtain instead of (8) the following equations from (4) and (14)

$`b_1(\omega )`$ $`=`$ $`a_1(\omega )`$
$`b_2(\omega )`$ $`=`$ $`a_2(\omega )+2\omega _\mathrm{p}c^1{\displaystyle \underset{k=\mathrm{}}{\overset{\mathrm{}}{}}}g_k(X_\mathrm{p}(\omega k\omega _\mathrm{q})+X_\mathrm{s}(\omega k\omega _\mathrm{q}))`$ (16)
Let suppose the local oscillator field in the form of

$$E_\mathrm{L}(t)=A_\mathrm{L}(t)\mathrm{cos}(\omega _\mathrm{p}t+\varphi )$$
(17)
where the dependence of the amplitude $`A_\mathrm{L}(t)`$ on $`t`$ is much slower than $`\mathrm{cos}\omega _\mathrm{p}`$t. Then for the envelope of the local oscillator field $`A_\mathrm{L}(t)`$ the Fourier expansion similar to (14) is valid

$$A_\mathrm{L}(t)=\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e_n\mathrm{exp}(\mathrm{i}n\omega _\mathrm{q}t)$$
(18)
The photodetector current has now the following form

$$I_{\mathrm{pd}}A_\mathrm{L}(t)(b_1(t)\mathrm{cos}\varphi +b_2(t)\mathrm{sin}\varphi )$$
(19)
and in the frequency domain one has

$$I_{\mathrm{pd}}(\omega )\mathrm{cos}\varphi \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e_nb_1(\omega n\omega _\mathrm{q})+\mathrm{sin}\varphi \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e_nb_2(\omega n\omega _\mathrm{q})$$
(20)
Let consider different parts in the photodetector output. The first term in equation (20) depends only on the amplitude fluctuations of the input field according to (16)

$$\mathrm{cos}\varphi \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e_nb_1(\omega n\omega _\mathrm{q})=\mathrm{cos}\varphi \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e_na_1(\omega n\omega _\mathrm{q})$$
(21)
The second term in equation (20) contains the signal and the noise parts. The noise part consists of the additive noise and the back action noise and has the following expression according to (15) and (16)

$`[\mathrm{sin}\varphi {\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e_nb_2(\omega n\omega _\mathrm{q})]_{\mathrm{noise}}=\mathrm{sin}\varphi {\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e_na_2(\omega n\omega _\mathrm{q})+2\omega _\mathrm{p}c^1\mathrm{sin}\varphi \lambda `$
$`{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}{\displaystyle \underset{k=\mathrm{}}{\overset{\mathrm{}}{}}}e_ng_kG(\omega k\omega _\mathrm{q}n\omega _\mathrm{q})\{{\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}g_ma_1(\omega k\omega _\mathrm{q}n\omega _\mathrm{q}m\omega _\mathrm{q})\}`$ (22)
Let consider only the photocurrent at small frequencies $`\omega \omega _\mu `$. Then the main input into the photocurrent will be given by the resonant terms for which $`k+n=0`$. With this supposition one has from equation (22)

$`[\mathrm{sin}\varphi {\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e_nb_2(\omega n\omega _q)]_{\mathrm{noise}}=\mathrm{sin}\varphi {\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e_na_2(\omega n\omega _\mathrm{q})+`$
$`\mathrm{sin}\varphi \lambda \xi (\omega ){\displaystyle \underset{m=\mathrm{}}{\overset{\mathrm{}}{}}}e_mg_m{\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}g_na_1(\omega n\omega _\mathrm{q})\}`$ (23)
Comparing equations (21) and (23) one can conclude that full compensation of back action noise is possible only for

$$e_n=\alpha g_n$$
(24)
where $`\alpha `$ is the same for all numbers $`n`$ so the forms of pump and local oscillator fields have to be the same (apart from the scale factor $`\alpha `$).

Let now consider the signal part of the second term in the r.h.s. of equation (20). From equations (7), (16) and (20) one has

$`[\mathrm{sin}\varphi {\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e_nb_2(\omega n\omega _\mathrm{q})]_{\mathrm{signal}}=`$
$`\mathrm{sin}\varphi {\displaystyle \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}}e_n\{{\displaystyle \underset{k=\mathrm{}}{\overset{\mathrm{}}{}}}g_k\xi (\omega k\omega _\mathrm{q}n\omega _\mathrm{q})F_\mathrm{s}(\omega k\omega _\mathrm{q}n\omega _\mathrm{q})\}`$ (25)
Evaluation of this expression for the condition $`k+n=0`$ gives

$$[\mathrm{sin}\varphi \underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e_nb_2(\omega n\omega _\mathrm{q})]_{\mathrm{signal}}=\mathrm{sin}\varphi \xi (\omega )F_\mathrm{s}(\omega )\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}e_ng_n$$
(26)
Combining equations (20), (23), (24) and (26) and supposing that the back action noise is compensated in the output of the photodetector one can obtain for the spectral density of noises in the photocurrent the following expression

$$N(\omega )\mathrm{sin}\varphi N_0\underset{n=\mathrm{}}{\overset{\mathrm{}}{}}g_ng_n=\mathrm{sin}\varphi N_0P$$
(27)
where it is supposed that fluctuations at frequencies $`\omega n\omega _q,n=0,1\mathrm{}`$ are uncorrelated and have the same spectral density $`N_0`$ (this assumption is valid for not very small duration of pump pulses), $`P`$ is proportional to the time averaged power of the pulsed pump. Then for the signal-to-noise ratio $`\mu `$ one has from equations (26) and (27) the following expression

$$\mu N_0^1P_{\mathrm{}}^{\mathrm{}}\xi (\omega )F_\mathrm{s}(\omega )^2d\omega $$
(28)
This value is just equal to the signal-to-noise ratio for continuous wave pump with a power $`P`$ and correlative processing of the output (cf. equation (11)). Note that the sensitivity here is not limited by the SQL like in the case of correlative processing of quadratures for the monochromatic pump and is increasing with the increase of $`P`$.

It is worth to mention that the condition for the back action noise compensation for the pulsed pump is just the same as for the monochromatic pump (cf. equation (12)) with substitution of the $`A^2`$ with the time averaged value $`P`$. Therefore the compensation of the back action noises for the finite frequency band can be possible for the time varying phase of the local oscillator .

## 6 Conclusion

The pumping of the displacement transducer with a train of the short high-intensity laser pulses is considered. The algorithm of optimal signal processing for such transducer is revealed. It consists of the correlative processing of the output using the pulsed local oscillator with the same envelope as for the pump field (apart from the scale factor). In this case the back action noise due to the radiation pressure force can be fully compensated and the sensitivity of the scheme to a detection of a classical external force is not limited by the SQL (as for the case of correlative quadrature processing and monochromatic pump field).

The pulsed pump can be advantageous over the single frequency pumping when the nonlinear optical elements are used unside the system. Thus considerable increase in sensitivity can be achieved for a gravitational interferometric Fabry-Perot type detector with a nonlinear optical element placed in a waist of the beam . The use of the phase-conjugate mirrors in a gravitational detector of the Michelson type allows to construct the system with the parallel arms . For such systems an efficiency depends on the instant power of the light beam and can be high for the short intensive pulses.

In this article only the problem of the force detection with known spectrum is considered. The reconstruction of unknown external force acting on the displacement transducer with the pulsed pump below the standard quantum limit will be considered elsewhere.