text:
# Two-phase flow in porous media: dynamical phase transition

## I Introduction

The field of two-phase flow in porous media is very rich on problems of a complex nature. These involve the study of properties on many different length scales and how to bridge the gap between the scales. Further, it is often possible to distinguish between the study of invasion processes and steady-state properties. The methods used in this field range from a set of experimental techniques, theoretical descriptions on different length scales to several numerical models or methodsDullien (1992); Sahimi (1995).

In the petroleum engineering tradition, the main source of experimental information about flow systems has for decades been displacement experiments on core samplesTiab and Donaldson (1996). However, for increased theoretical or fundamental understanding, well controlled laboratory experiments have been the key. In particular, in the eighties, the experiments and simulations by Lenormand *et al.* lead to the development of a phase diagram for drainage, i.e., the displacement of a wetting phase by a nonwetting phaseLenormand et al. (1983, 1988). Drainage was classified into the regimes stable displacement, viscous fingering, and capillary fingering. The regime boundaries depend on two dimensionless numbers, the capillary number Ca and the viscosity ratio $`M`$.

Two-phase flow consists of more than pure displacement processes. The complexity of ganglion dynamics and the existence of different regimes of flow of bubbles and blobs has been investigated experimentally by the Payatakes groupAvraam and Payatakes (1995a, b, 1999). They used etched glass networks for their studies, varying flow parameters within large ranges. By simultaneous injection of two fluids the experiments examine steady-state properties or close to steady-state properties.

We use a simulator for two-phase flow in porous media that is based on the Washburn equationWashburn (1921), see Sec.II. This line of modeling dates back to the mid eighties to Koplik *et al.*Koplik and Lasseter (1985). The work by the Payatakes group is a continuation of this traditionConstantinides and Payatakes (1991, 1996). Most of the work on two-phase flow simulation has been in the actual two-phase regime. That is to say that pure invasion processes like drainage were simulated until breakthrough of the nonwetting phase. The aforementioned regimes, the phase diagram, of drainage are well established and their properties much studied. When it comes to simulations for investigating steady-state properties, the majority of the work involves finding relative permeability curves for a wide range of parameters. Much effort has been put into making simulators that are specific for a given porous medium and fluid system. Further, other simulation techniques exist. On a even more detailed scale than our modeling is the work using lattice Boltzmann methodsRothman (1990); Ferreol and Rothman (1995); Langaas and Papatzacos (2001).

We wish here to take one step back and take a broader look at these systems. Assuming a porous medium with two phases, one wetting and one nonwetting with respect to the medium, there are three possible states of flow: only the nonwetting phase flows, only the wetting phase flows, or both flow simultaneously. Depending on a large set of parameters the system will find itself within one of these situations. Changing system parameters within appropriate ranges causes the system to undergo what is called a dynamical phase transition. Using the language of thermodynamics and critical phenomena we provide in this paper a phase diagram for steady-state flow. As was the case for the phase diagram of Lenormand *et al.* the capillary number and the viscosity ratio are the parameters of interest also in this primary study of the steady-state phase diagram.

We use a not too complex network simulator that can do real steady-state simulations. Nevertheless, the resulting phase diagram is complex. We believe that the structure of this diagram is of general interest. Phase boundaries are located by the simulations. Further, based on the simulations a semi-empirical theory for the location of the phase boundaries is given. We believe the overall structure of the phase diagram is universal for the two-phase flow system, whereas the exact positions of phase boundaries will vary from system to system.

## II Model

In order to investigate question of interest regarding two-phase flow in porous media it is possible to apply methods that count as theoretical, experimental or numerical. Experimental methods provide quantitative ‘correct’ results for the actual system investigated. However, results may differ from sample to sample or from system to system. By means of simulations the system properties can be controlled to a greater extent. This is highly advantageous when one wishes to study the effect of varying parameters of the system.

The results of the paper are based on a network simulator for immiscible two-phase flow. This line of modeling which is based on the Washburn equationWashburn (1921) dates back to the work of several groups in the mid 1980’sKoplik and Lasseter (1985); Constantinides and Payatakes (1991, 1996). Our model is a continuation of the model developed by Aker *et al.*Aker et al. (1998a, b). Although a thorough presentation can be found in Knudsen et al. (2002), we provide for clarity a brief résumé of the main aspects.

The porous media are represented by networks of tubes, forming square lattices in two dimensions (2D) and cubic lattices in 3D, both tilted $`45^{}`$with respect to the imposed pressure gradient and thus to the overall direction of flow. We refer to the the lattice sites where four(2D) and six(3D) tubes meet as nodes. Volume in the model is contained in the tubes and not in the nodes, although effective node volumes are used in the modeling of transport through the nodes. For further details, seeKnudsen et al. (2002). Randomness is incorporated by distorting the nodes randomly within a circle (2D) and a sphere (3D) around their respective lattice positions. This gives a distribution of tube lengths in the system. Further, the radii are drawn from a flat distribution so that the radius of a given tube is $`r(0.1l,0.4l)`$, where $`l`$is the length of that tube.

The model is filled with two fluid phases that flow within the system of tubes. The flow in each tube obeys the Washburn equationWashburn (1921), $`q=(\sigma k/\mu )(\mathrm{\Delta }pp_c)/l`$. With respect to momentum transfer these tubes are cylindrical with cross-section area $`\sigma `$, length $`l`$, and radius $`r`$. The permeability is $`k=r^2/8`$which is known for Hagen-Poiseuille flowGerhart et al. (1992). Further, $`\mu `$is the viscosity of the phase present in the tube. If both phases are present the volume average of their viscosities is used. The volumetric flow rate is denoted by $`q`$and the pressure difference between the ends of the tube by $`\mathrm{\Delta }p`$. The summation is the sum over all capillary pressures $`p_c`$within the tube. With respect to capillary pressure the tubes are hour-glass shaped, meaning that a meniscus at position $`x(0,l)`$in the tube has capillary pressure $`p_c=(2\gamma /r)[1\mathrm{cos}(2\pi x/l)]`$, where $`\gamma `$is the interfacial tension between the two phases. This is a modified version of the Young-Laplace lawDullien (1992); Aker et al. (1998a).

Biperiodic (triperiodic in case of 3D) boundary conditions are used so that the flow is by construction steady flow. That is to say, the systems are closed so both phases retain their initial volume fractions, i.e., their saturations. These boundary conditions are in 2D equivalent to a flow restrained to be on the surface of a torus. The flow is driven by a globally applied pressure gradient. Usually invasion processes are driven by setting up a pressure difference between two borders, inlet and outlet. Since, by construction, the outlet is directly joined with the inlet in our system, we give instead a so-called global pressure drop when passing a line or cut through the systemRoux (unpublished); Knudsen et al. (2002). This is equivalent to imposing constraints on the pressure gradient that is experienced throughout the network. Integration of the pressure gradient along an arbitrary closed path one lap around the system, making sure to pass the ‘inlet-outlet’-cut once, should add up to the same global pressure difference.

The system is forward integrated in time using the Euler scheme. For each time step the distribution of the phases leads to a recalculation of the effective viscosities in the tubes and the capillary pressures across the menisci, and thus the coefficients in the equations for the pressure field. When the pressure field is known, the flow field follows automatically and the integration step may be performed. Details are found inKnudsen et al. (2002).

## III Phase diagrams

We present here phase diagrams for the dynamical properties of steady two-phase flow. The structure of the phase diagram is richer than one might naively believe. First we introduce the well-known concept of fractional flow (Eq.(2)) and provide a basic description of the nature of the two-phase flow system. We determine phase boundaries in the parameter space of saturation $`S_{\mathrm{nw}}`$ (defined in Eq.(1)), capillary number Ca (defined in Eq.(4)), and viscosity ratio $`M`$(defined in Eq.(3)). The definition of Ca is a crucial point since the definition is not unique. For the numerical exploration of the phase diagram one definition is used consistently in all simulation series: first at constant Ca and varying $`S_{\mathrm{nw}}`$, second at constant $`S_{\mathrm{nw}}`$and varying Ca. While Sec. III contains simulated diagrams, Sec. IV readdresses the definition of the capillary number and contains a semi-empirical theory for the location of the phase boundaries that are obtained in Sec. III. Upon presenting the method of finding the phase boundaries, the order of the transition is briefly discussed. However, the final discussion of the order is postponed to Sec. V.

### III.1 Order parameter

In steady flow the volume fractions of the wetting and the nonwetting phase will not change in time. In our simulation model this requirement is fulfilled by having a closed system. The nonwetting saturation of the system is defined as

$$S_{\mathrm{nw}}=\frac{V_{\mathrm{nw}}}{V_{\mathrm{tot}}},$$
(1)
where $`V_{\mathrm{nw}}`$and $`V_{\mathrm{tot}}`$are the nonwetting and total volume, respectively. The wetting saturation $`S_\mathrm{w}`$is defined similarly. For each value of the saturation one can in the simulations fix the total flux $`Q_{\mathrm{tot}}`$in the system. In order to keep the total flux constant the imposed pressure gradient fluctuates in time around some typical mean value. This mean pressure is a function of saturation and total flux. Further, one can measure the flux of each of the phases in the system, $`Q_{\mathrm{nw}}`$and $`Q_\mathrm{w}`$, respectively. One defines the nonwetting fractional flow as

$$F_{\mathrm{nw}}=\frac{Q_{\mathrm{nw}}}{Q_{\mathrm{tot}}},$$
(2)
and likewise for the wetting fractional flow. A sample of nonwetting fractional flow as a function of nonwetting saturation is given in Fig. 1. The transport properties of the system depend on dimensionless groups of system parameters, namely the viscosity ratio and capillary number. We define the viscosity ratio as the nonwetting viscosity divided by the wetting viscosity;

$$M=\frac{\mu _{\mathrm{nw}}}{\mu _\mathrm{w}}.$$
(3)
Further, the capillary number is defined as

$$\mathrm{Ca}=\mathrm{Ca}_{\mathrm{dyn}}=\frac{Q_{\mathrm{tot}}\mu }{\gamma \mathrm{\Sigma }},$$
(4)
where $`\gamma `$is the interfacial tension between the two phases, $`\mathrm{\Sigma }`$ is the cross-section area of the system, and $`\mu `$is the effective (weighted) viscosity of the system; $`\mu =S_{\mathrm{nw}}\mu _{\mathrm{nw}}+S_\mathrm{w}\mu _\mathrm{w}`$. Physically the capillary number measures the ratio between viscous and capillary forces in the system. Using the weighted viscosity is motivated by the fact that when both phases flow, both viscosities play a role in determining the relative strength of viscous and capillary forces. We use this *dynamical* definition of Ca in this section, upon the exploration of the phase space. For convenience we leave out the subscript of $`\mathrm{Ca}_{\mathrm{dyn}}`$, but we will use it in the subsequent discussion in Sec. IV for clarity.

As illustrated in Fig. 1 fractional flow curves depend on $`M`$and Ca. This dependence is nontrivial. The curves exhibit three different regimes. For low nonwetting saturation only the wetting phase flows. Likewise, for high nonwetting saturation only the nonwetting phase flows. These two regions are thus effectively single-phase flow in a constraining environment consisting of both the solid porous medium and the immobilized phase. In-between these regions is the two-phase flow region. The cross-overs from single-phase flow to two-phase flow are in fact dynamical phase transitions. We wish to study this system with that perspective. As order parameter for the transition from single wetting flow to two-phase flow we use the nonwetting fractional flow. Similarly we take the wetting fractional flow as the order parameter for the transition from single nonwetting flow to two-phase flow. In this way the order parameter is zero for single-phase flow and nonzero for two-phase flow.

Having defined the order parameter, the nature of possible dynamical phase transitions in the system is well illustrated in Fig. 1. Starting at the left-hand-side, the transition is from single-phase wetting flow to two-phase flow when increasing the nonwetting saturation. Both samples show a continuous transition. That is to say that at least to the resolution of the curves the order parameter changes continuously from zero to nonzero values at the phase boundary. The transition of the curve marked with $``$on the right-hand-side has clear signs of being first order. Not only can we see how the order parameter is discontinuous, but the one isolated data point of $`F_\mathrm{w}=0.0`$indicates the presence of hysteresis or history dependence in the system. Hysteresis is typical for first order transitions. In fact, we will use studies of hysteresis in order to unravel the nature of the phase diagram. As to the other curve, marked with $`\times `$, on the right-hand-side, the transition seems to be continuous or weakly first order.

### III.2 Hysteresis loops

In order to map out the phase diagram, we need a tool that can locate the phase boundaries in parameter space. The order of the transitions should also be determined. Simulations of single points in parameter space to steady-state as in Fig. 1 are good for the part of parameter space where history effects are small. However, whenever one is close to a first order transition the history may give data points indicating single-phase flow or two-phase flow. It is possible to run simulations on several statistically equal realizations of the porous media, as well as from various initial configurations, in order to find the region where both single-phase and two-phase flow are possible. The computations are somewhat computer demanding, and we have found it more efficient to simulate entire hysteresis loops to this effect. For instance starting at a specified saturation at which the system has two-phase flow, we change in a stepwise fashion the saturation slightly, run until steady-state, and measure the fractional flow. This is done until the system passes the transition to single-phase flow. Thereafter the saturation is changed the opposite way until the cross-over from single-phase flow to two-phase flow occurs. In case of a continuous transition, these two points of transition will be the same. However, in case of a first order transition these points will not be the same, and the fractional flow will appear as a hysteresis loop, see Fig. 2(a).

It is also possible to keep the saturation constant and vary some other parameter, typically the capillary number. Similar loops are expected at first order transitions by this procedure. The advantage of changing Ca is that one simply has to change the total flux $`Q_{\mathrm{tot}}`$in the system. This is a unique operation, while changing the saturation smoothly is not unique. It requires that one makes a choice as to where in the system the changes should be made. In general, we increase or decrease existing bubbles by moving the menisci defining their surface. If the system is very fragmented, some smaller bubbles are typically also removed in order to reach the desired new saturation.

### III.3 Diagrams for constant Ca

In this subsection we study one sample of the porous network of size $`20\times 40`$. We do not yet go into the possible dependence on size and topology, but focus on how the results depend on the three parameters: phase saturation, capillary number Ca, and viscosity ratio $`M`$. First we keep Ca and $`M`$constant while we vary the saturation stepwise over the transition from single-phase flow to two-phase flow and back.

Examples of transitions between single-phase and two-phase flow are shown in Fig. 2. By varying the saturation both ways over the transition, we observe how there can (a), and cannot (b), be hysteresis in the transition. Having a transition with little or without hysteresis as in Fig. 2(b), leads to a quite precise determination of the transition point. On the other hand it is not so clear where the transition occurs for the sample in Fig. 2(a). Not only are there two actual transition points, but these two points should be expected to depend on the actual path in parameter space. By that we mean in this case how many steps there are in saturation, and for what time the system is allowed to relax at each step. Further, the variation from sample to sample, or from realization to realization, of systems that are in principle similar, may very well be larger when there is hysteresis of this kind. For curves as in Fig. 2(b), the statistical variations are smaller. The transitions in Fig. 2 were taken to be at (a) $`S_{\mathrm{nw}}=0.26`$ and (b) $`S_{\mathrm{nw}}=0.84`$. Here we employed that convention that the system is single-phase whenever it is single-phase in one of the simulated directions. This definition is mainly motivated by the need of consistency in the extraction of data points in the following. This is in contrast to the more conventional definition, where the transition is located at the steepest point of the order parameter curve. However, this definition leeds to too much noise in the analysis.

The range of capillary numbers studied is $`\mathrm{Ca}(3.16\times 10^2,3.16\times 10^4)`$. On a base 10 logarithmic scale 9 values of Ca are chosen in steps of $`0.25`$from $`1.50`$ to $`3.50`$. Similarly 21 values of $`M`$are chosen in steps of $`0.25`$in the interval $`\mathrm{log}_{10}M(2.50,2.50)`$, which means that $`M(3.16\times 10^3,3.16\times 10^2)`$.

The resulting phase boundaries are shown in Fig. 3 as data points. The solid lines are constructed phase boundaries following from the semi-empirical theory given in Sec. IV. For now, let us consider the lines as visual aids to better see the phase boundaries in the diagrams. With increasing nonwetting saturation, the two lines separate the dynamical phases: single-phase wetting flow, two-phase flow, and single-phase nonwetting flow. Observe that for high capillary numbers the two-phase region spans almost all of the saturation range. As the capillary number is lowered, the single-phase regions grow in size. There is some kind of symmetry between the two transitions. That is to say the same effects occur on both sides, although not at the same numerical values.

In order to treat the two transitions simultaneously, we employ the following notions. Since at each transition there is more volume of the phase that is close to single-phase flow, we refer to this phase as the majority phase. The other phase is the minority phase. On the left side the majority phase is wetting and the minority phase is nonwetting. Regarding viscosity ratios $`M`$there are basically three possibilities: viscosity matching phases, favorable viscosity ratio, and unfavorable viscosity ratio. Considering the system from the point of view of the majority phase, the viscosity ratio is favorable when the majority phase is the less viscous phase. Likewise, the viscosity ratio is unfavorable when the minority phase is less viscous. On the $`y`$-axes of Fig. 3 $`M`$varies over more than four decades. The positive side (on the logarithmic scale) is when the nonwetting phase is more viscous than the wetting phase. Thus, the viscosity ratio is favorable on the left side and unfavorable on the right side. For negative values of $`\mathrm{log}_{10}(M)`$the ratio is unfavorable on the left side and favorable on the right side.

The general structure of the phase diagrams in Fig. 3 is as follows. First, consider the case of viscosity matching phases, $`M=1`$. Here, we observe that for high capillary numbers the two-phase region spans almost all of the saturation range. As Ca is lowered the single-phase regions on each side open up, and the two-phase region becomes smaller. We cannot tell from the presented diagrams the limiting behavior of the phase boundaries for small Ca. However, since the nine diagrams in Fig. 3 are evenly distributed on the log-scale, one can conclude that, roughly speaking, the dependence on Ca is logarithmic(also cfr. Fig. 4). It is reasonable to expect this behavior to persist also for lower Ca.

The dependence on $`M`$is interesting. For all values of Ca the fact that each curve has two main parts persists. Roughly, there is one situation for favorable viscosity ratios and another for unfavorable viscosity ratios. The subsequent discussion of these two situations are made from the two-phase flow point of view. The single-phase point of view is taken in Sec. IV.

For unfavorable viscosity ratio (lower left and upper right sides in each of the diagrams in Fig. 3) the phase boundaries are almost constant with respect to the value of $`M`$. As long as there are clear majority and minority phases, and one phase is close to going single-phase, the minority phase consists of bubbles and ganglions that are held back by capillary forces. It is always the case that the minority phase is held back in this way. However, when the viscosity ratio is unfavorable, the minority phase flows more easily. The viscous drag on the minority phase is thus negligible for this range of $`M`$-values explaining the constant location of the boundaries. This can be understood so that when the minority phase flows sufficiently easy, then it simply follows the majority phase around when capillary forces are overcome. Notice that the noise level increases with decreasing Ca. When the capillary forces become relatively stronger the motion of bubbles involves larger fluctuations, so this behavior is expected when properties are averaged over the same amount of time.

A different situation appears for favorable viscosity ratios. The capillary forces and viscous forces both try to hold back the minority phase. In this situation the actual viscosity ratio plays a role. The more favorable viscosity ratio, the more will the minority phase slow down. Because of the capillary forces the minority phase will not only slow down, but actually stop, and single-phase flow results. The border between two-phase flow and single-phase flow varies with $`M`$, as can be seen in Fig. 3. We observe the general trend that the single-phase region increases logarithmically with increasingly favorable viscosity ratio.

### III.4 Diagrams for constant $`M`$

So far we have discussed the results from simulations where the saturation has been varied over the dynamical phase transition. In this subsection we present results from simulations in which the saturation and the viscosity ratio are held constant, while Ca is varied down and up again over the transition. Because of history effects and the fact that saturation and Ca are not varied in the same way, it is a priori not clear that the phase boundaries will be the same. Therefore, we check the results for consistency.

For nine different values of $`M`$, equidistantly distributed on a logarithmic scale between $`M=1/100`$to $`M=100`$, simulations have been performed for a large set of saturations. The results are shown in Fig. 4. For large viscosity contrasts it turned out to be difficult to simulate very small saturations of the minority phase, possibly because the transition would occur for very large capillary numbers. More importantly, the central part of the saturation range was also difficult to simulate, because the fluctuations in the systems are larger in this region. One might expect that it is possible for the system to reach either single-phase wetting flow or single-phase nonwetting flow simply as the result of fluctuating into that state. We found that the fluctuations remain large for capillary numbers all the way down to $`\mathrm{Ca}=10^6`$for several samples. It is also possible that there is a finite saturation range, where the dynamical phase is two-phase flow for all values of Ca. We conclude that this region therefore may be out of reach. Even if we could get the results from the simulator, such low capillary numbers are outside the range where this simulator works well, and the relevance to physical systems would be limited.

The observations we can make from Fig. 4 are as follows. The nine presented diagrams have basically the same appearance. They consist of three parts: single-phase wetting flow in the lower left corner, single-phase nonwetting flow in the lower right corner, and two-phase flow in the upper middle region.

First, observe that as far as simulations have been performed, the phase boundaries have ending points at zero and unity saturation for finite values of Ca. At higher values of Ca the system has two-phase flow for all saturations. This is not surprising. High capillary numbers mean vanishing small capillary forces relative to viscous forces. In principle, as long as the viscosity contrast is finite or not too large, the more viscous phase would still flow unless there also are capillary forces to keep it from flowing. Hence, for these viscosity ratios, the capillary forces are too small above some threshold to prevent one phase from flowing.

Second, the general behavior of the two boundaries away from their end points, is that as Ca is lowered, the minority saturation increases. The meaning is simply that when the capillary forces become relatively stronger the single-phase regions grow in size. For low Ca values the two-phase region becomes quite narrow. As we already have pointed out, we cannot tell whether the two boundaries meet in a triple point or if there will be a finite two-phase region also for smaller Ca.

In order to compare the curves in Fig. 4, the results have been extracted and put into one diagram, see Fig. 5(b). The curves overlap considerably so we have chosen to use the same line style for most of the lines. The exception is that for the three most favorable viscosity ratios of each transition a different line style is used. The three curves separate more and more from the rest of the curves with increasing viscosity contrast.

To compare the simulation results shown in Figs. 3 and 4, we exhibit data from Fig. 3 in Fig. 5(a). Here, the same set of viscosity ratios were used as in Fig. 5(b). This was done for best possible comparison. The line styles are also the same. Note that the range of capillary numbers in Fig. 5(a) is smaller than in Fig. 5(b). The results are not identical, but nevertheless quite close to each other within the overlapping range. We conclude that the two ways of performing the simulations give consistent results.

There are some general features of the diagrams in Fig. 5. First we observe that for both the left side boundary and the right side boundary, there seems to be a lower limit. This means that for all viscosity ratios, the respective single-phase region on each side extends at least to this limit. Roughly speaking these limits vary logarithmically with Ca. Again we cannot tell if this behavior will continue for lower Ca than depicted.

The other feature is that for higher Ca, say $`\mathrm{Ca}>10^{3.5}`$, favorable viscosity ratio affects the phase boundary. To be precise, for the three intermediate values of $`M`$: $`M=10^{0.5}`$, $`M=1`$, and $`M=10^{0.5}`$, the phase boundaries on both sides are more or less equal to the limiting behavior just discussed. For these values of $`M`$, $`M=10^{1.0}`$, $`M=10^{1.5}`$, and $`M=10^{2.0}`$, the viscosity ratio is favorable for the transition on the left side between single-phase wetting flow and two-phase flow. The corresponding three lines are marked with dashed lines. These values of $`M`$are unfavorable for the right-hand side transition. The last three values of $`M`$, $`M=10^{1.0}`$, $`M=10^{1.5}`$, and $`M=10^{2.0}`$, are favorable for the right-hand side transitions. They are indicated with dash-dotted lines. For the left-hand side these $`M`$values are unfavorable. In conclusion, there are six lines that are singled out for their favorable viscosity ratio. They clearly show how the single-phase regions become larger as the viscosity ratio becomes increasingly favorable. This effect can be seen in both Figs. 5(a) and 5(b).

## IV Location of the phase boundaries

In the previous subsection the point of view taken was that of the two-phase flow system. That is to say, we used the dynamical definition of the capillary number, Eq. (4), and we asked when does one of the phases cease to flow. Reversing the point of view is fruitful. Consider the system starting from a single-phase flow situation, the natural question is when does the second phase start to flow.

Being in a single-phase flow situation, although with both phases present in the system, means that the capillary number should be defined differently. Namely, the relevant viscosity is no longer the volume averaged viscosity, but the viscosity of the fluid that is actually flowing. The wetting and nonwetting capillary number is defined as

$$\mathrm{Ca}_\mathrm{w}=\frac{Q_{\mathrm{tot}}\mu _\mathrm{w}}{\mathrm{\Sigma }\gamma },$$
(5)
and

$$\mathrm{Ca}_{\mathrm{nw}}=\frac{Q_{\mathrm{tot}}\mu _{\mathrm{nw}}}{\mathrm{\Sigma }\gamma },$$
(6)
respectively. Here, $`\mathrm{Ca}_\mathrm{w}`$ is the relevant capillary number for single-phase wetting flow, whereas, $`\mathrm{Ca}_{\mathrm{nw}}`$ is relevant for single-phase nonwetting flow.

As opposed to the qualitative statements regarding the phase boundaries that were made in Sec. III, these definitions of the capillary number are also quantitatively useful. For a given system, i.e. given topology and saturation, the transition to two-phase flow takes place at a given wetting (or nonwetting) capillary number. That means with varying value of the viscosity ratio $`M`$, which also changes $`\mathrm{Ca}_{\mathrm{dyn}}`$, the relevant single-phase capillary number still holds a fixed value at the transition. This is plausible as long as the transition to two-phase flow is continuous. In that case the amount of minority phase flowing just after its mobilization is still only a small fraction of the total flow, and therefore the physical effective viscosity is still very close to the viscosity of the single-phase that just before flew alone. However, should the transition be first order so that a finite fraction is mobilized at once, then for large variations in $`M`$ some corrections should be made. Nevertheless, for now we assume that the single-phase capillary number governs the location of the phase boundaries. The question of the order of the transition is addressed in Sec. V, and thereby the validity of this assumption is discussed.

We can write the dynamical capillary number

$$\mathrm{Ca}_{\mathrm{dyn}}=\frac{Q_{\mathrm{tot}}\mu _\mathrm{w}}{\mathrm{\Sigma }\gamma }\left(S_\mathrm{w}+MS_{\mathrm{nw}}\right),$$
(7)
from which it follows that in the limit of small $`M`$

$$\mathrm{Ca}_{\mathrm{dyn}}=\mathrm{Ca}_\mathrm{w}S_\mathrm{w},$$
(8)
and in the limit of large $`M`$

$$\mathrm{Ca}_{\mathrm{dyn}}=\mathrm{Ca}_{\mathrm{nw}}S_{\mathrm{nw}}.$$
(9)
In Fig. 5 we observe the limiting behavior for unfavorable $`M`$ of the phase boundaries. There is on each side a limiting line, below which there is always single-phase flow. By Eqs. (8) and (9) these limits are identified with the certain wetting or nonwetting capillary number that for each value of the saturation determines the location of the transition. The limiting lines and the fitted functional relationships are shown in Fig. 6. We take this logarithmic dependency as an empirical input, and see which predictions can be made therefrom employing theory.

We want to calculate the location of the phase boundaries for fixed dynamical capillary numbers as shown in Fig. 3. Considering first the left hand side (transition from wetting single-phase to two-phase flow), the criterion for the phase boundary is expressed as

$$\mathrm{Ca}_\mathrm{w}S_\mathrm{w}=\mathrm{Ca}_l10^{\alpha _lS_{\mathrm{nw}}},$$
(10)
thus being a function of saturation. Letting $`\mathrm{Ca}_{\mathrm{dyn}}`$ be fixed, this criterion inserted into Eq. (7) gives when solving for $`M`$ as a function of $`S_{\mathrm{nw}}`$

$$M=\frac{1S_{\mathrm{nw}}}{S_{\mathrm{nw}}}\left(\frac{\mathrm{Ca}_{\mathrm{dyn}}10^{\alpha _lS_{\mathrm{nw}}}}{\mathrm{Ca}_l}1\right).$$
(11)
The interpretation of Eq. (11) is that for a fixed value of $`\mathrm{Ca}_{\mathrm{dyn}}`$ the phase boundary is given by the set of values of $`M`$ and $`S_{\mathrm{nw}}`$ that fulfill this relationship. For values of $`\mathrm{Ca}_{\mathrm{dyn}}>\mathrm{Ca}_l`$, the wetting single-phase region does not exist, and so does the identification in the limit, Eq. (8) not apply. Still the criterion Eq. (10) can be used, and the result is that the single-phase region can only be reached for values of $`M`$ larger than a certain threshold, see Fig. 3(a)-(c). This is not the case for $`\mathrm{Ca}_{\mathrm{dyn}}<\mathrm{Ca}_l`$, where there is a single-phase region for all $`M`$.

In a similar way the boundaries on the right hand side between two-phase flow and single-phase nonwetting flow is found to obey

$$M=\left[\frac{1S_\mathrm{w}}{S_\mathrm{w}}\left(\frac{\mathrm{Ca}_{\mathrm{dyn}}10^{\alpha _rS_\mathrm{w}}}{\mathrm{Ca}_r}1\right)\right]^1.$$
(12)
By using the Eqs. (11) and (12) the solid lines in Fig. 3 are the constructed phase boundaries for the respective values of $`\mathrm{Ca}_{\mathrm{dyn}}`$. The general good agreement between the simulated data points and these semi-empirical, semi-theoretical lines, serves as a confirmation a posteriori of this theoretical approach. Indeed, we conclude that the single-phase capillary number is the relevant entity that controls the transition to two-phase flow, rather than the more commonly used capillary number defined in Eq. (4). However, we should point out that from an experimental point of view, it is the conventionally defined capillary number which is the most relevant one.

## V Order of the transitions

So far we have only discussed the location of the dynamical phase boundaries. The order of the transitions across these boundaries should also be discussed. All the simulations presented so far were done on a single system size, namely $`20\times 40`$nodes. It is necessary to take into account effects of changing the system size to determine the order of the transitions. However, first we present results for this system size.

Each simulation in Fig. 3 can be categorized into showing or not showing hysteresis. We look for the general trend while accepting that fluctuations and randomness play a role. Some realizations do not show hysteresis even though all neighboring points do, and vice versa. It turns out that the phase boundaries of Fig. 3 can be divided into two logical parts. The almost vertical line segments where the systems have unfavorable viscosity ratio differ from the other parts of the curves where the viscosity ratio is favorable. In the vertical parts the simulations show no or little signs of hysteresis in the order parameter. This indicates that the transitions in these regions are continuous. The fact that the curves are noisy for lower Ca in these regions, does not affect the order. In the sloping parts of the boundaries the order parameter exhibit discontinuity and hysteresis in most simulations. The level of hysteresis increase with increasing viscosity contrast. We claim that the transitions in these regions are first order.

Whether a transition is continuous or first order is a property connected to a specific path in parameter space. However, the paths in Fig. 3 and in Fig. 4 all mean a transition from single-phase flow to two-phase flow. We find similar hysteresis effect for both kinds of paths through parameter space.

In general, it is possible that the appearance of different orders of the transitions is a result of finite size effects. We have performed some selected simulations on a larger system size; $`40\times 80`$nodes. Comparisons between the results of the two system sizes are found in Fig. 7. The first order transition in (b) persist also for the larger system. If the signs of first order had been weaker for larger system sizes, then finite size effects had been the cause for apparent first order transitions. However, this is not the case here.

The noisy character of the vertical parts, showing continuous transitions, can be attributed to the geometrical heterogeneity of the system. Some parts of the system are more active in transportation than others. In such rather small systems self-averaging is small and by chance the minority phase is to a larger or smaller extent placed in inactive parts of the system. This gives a shift in position of the transition, i.e. noise, but does not affect the correlation between minority phase bubbles that are placed in active regions of the system. In fact these bubbles that are only held back by capillary forces are all correlated upon onset of mobilization. The motion of one bubble in one place relieves the pressure on all other minority bubbles in the system. The opposite is also true, when this bubble stops again, the pressure on the others increase immediately and this is likely to cause the onset of mobilization elsewhere. This is consistent with having infinite correlation length in statistical mechanics at continuous transitions.

The sloped sections of the boundaries, which from the hysteresis loop consideration show signs of first order transitions, show different characteristics. The viscosity ratio is unfavorable and the mobilization requires more pressure on a bubble, but when mobilized, larger bubbles or connected bubbles are involved. However, this effect is more local in the sense that the minority phase is not moving a little bit here and a little bit there in the form of different bubbles, but rather one larger region moves for a longer time when it is first mobilized. This correspond to a non-infinite correlation length. Note that the system is driven by the global criterion of fixed total flux through the system. From this point of view, everything is all the time correlated with everything, so when we speak about correlations, we speak about the motions of the minority phase as they are observed.

## VI Conclusion

This study concerns steady-state properties of two-phase flow in porous media. Average flow properties are monitored by using a network simulator on pore level. A systematic change of the system parameters: phase saturation, viscosity ratio, and capillary number, is performed. We demonstrate how the system can be in either one of three dynamical states or phases. These are: single-phase wetting flow, two-phase flow, and single-phase nonwetting flow. Upon passing from one dynamical phase to another the system undergoes a dynamical phase transition. The phase diagram for the dynamical phases is revealed and its properties are discussed.

By means of hysteresis we determine the order of the phase transitions and we find that both continuous and first order transitions occur across different parts of the phase-boundaries. This is connected to two very different situations. One situation is when the viscosity ratio is defined to be favorable with respect to the majority phase: both capillary forces and viscous forces hold back the minority phase from flowing, next to this transition. The other case is when the viscosity ratio is defined to be unfavorable: only capillary forces try to hold back the minority phase next to the transition.

The phase diagram that is provided is for steady-state properties. This should be set in perspective to the phase diagrams by Lenormand *et al.* for drainage invasion flow propertiesLenormand et al. (1983, 1988). Although the fact that the system investigated is a particular steady-state two-phase flow system, we argue that the general structure of the dynamical phase diagram has a universally valid structure.

We find that the definition of the capillary number is of major importance. Whereas for the actual two-phase flow region we employ a volume averaged viscosity in the definition, it turns out that the actual locations of the phase boundaries are determined by the capillary number based on the viscosity of the single phase in question. Based on the relationships between the capillary numbers, theoretical part, and limiting of the phase boundaries, found by the simulations, we establish a semi-empirical theory for the location of the phase boundaries. Upon direct comparison this theory is found to be in agreement with the simulated data points.

###### Acknowledgements.

H.A.K. thanks VISTA, a collaboration between Statoil and the Norwegian Academy of Science and Letters, for financial support. We thank E. Skjetne for very valuable discussions, in particular in connection with Section IV where his contribution was substantial.


format: text

text:
# Three dimensional imaging of short pulses

## 1 Introduction

The study of ultrafast phenomena has been a major scientific priority during last decades covering different topics such as the study of radiation-matter interactions hullier , transient response of molecules and atoms zewail , coherent control of chemical reactions Baumert  or communication and information technology stegeman . The growth of this field relies upon the development of sources of femtosecond radiation and of appropriate techniques able to provide time domain information in the femtosecond scale. However, during the interaction of short optical pulses with a nonlinear medium, different mechanisms can lead to its reshaping into complex spatio-temporal structures with non-trivial light distribution PDT00 . As a consequence, their complete characterization requires a method capable of acquiring a snap-shot of their intensity distribution in the whole 3-dimensional (3D; $`x,y,tz/c`$) comoving frame.

Most of the available methods for pulse diagnostic provide information of the WP characteristics in a space of reduced dimensionality. The use of frequency resolved autocorrelation techniques (i.e. FROG, SPIDER) allows for the recovery of the temporal intensity and phase profile of a given pulse but assumes uniform transverse spatial distribution Trebino ; Jaconis . On the contrary, the characterization of transversally localized beams often relies upon the optical imaging onto CCD cameras, therefore the temporal information is lost because of their integration times unavoidably larger than the optical pulse duration. Recently, a space-time characterization method based on extended SPIDER technique has been developed capable of resolving electric field characteristics in time and along one spatial coordinate dorrer . A quite direct way of obtaining spatio-temporal intensity profiles of a WP is to perform measurements with a streak camera, which allows a temporal resolution up to fractions of ps streak . This technique allowed the investigation of the dynamics of the breakup along the pulse envelope of a large elliptical beam propagating into a saturable Kerr nonlinear medium Lantz02a ; Lantz02b . However, also in this case, the space-time maps are intrinsically two dimensional (1 spatial + temporal dimensions).

A different approach to the problem considers the retrieval of the pulse shape through an all-optical processing by means of spatially resolved detection systems combined with gating techniques. The principle of the method is that of characterizing with spatial resolution an optical field that is proportional to the product $`E_O(𝐱,t)E_R(𝐱,t)`$, where $`E_O(𝐱,t)`$ is the object to be measured and $`E_R(𝐱,t)`$ is a suitable reference pulse. Since the product is different from zero only on the intersection of the support of both fields111Rigorously, the support of a function is the set of points where its value is different from zero. For realistic optical fields that have exponentially decaying tails we can define a “practical” support, defined as the set of points in the space-time in which the field amplitude is larger than $`1/e`$ times the peak value., by translating the reference with respect to the object, we get the possibility of recording information from different parts of the object. Among the linear time gating techniques, light–in–flight holographic recording has been the first technique which permitted the recording of dynamically evolving light fields during propagation Denisyuk69 ; Abramson83 ; Abramson89 . Recently, this technique has been adapted to study the propagation of a 3 ps long pulse in linear media Kubota02 . Linear probing techniques were also exploited in order to obtain time resolved imaging, like the probing of the birefringence properties of plasma by means of a delayed, spatially-extended 100 fs pulses to investigate the dynamics of laser pulse focusing in air Fujimoto99 .

Nonlinear processes have been employed since long ago to resolve in time the evolution of ultrafast phenomena. Among them, the quadratic nonlinearity has been proved to be particularly versatile due to the fact that it provides easily terms containing the product of two optical fields. Recently, a type II degenerate parametric amplification scheme has been employed to obtain time resolved 2D images of a ps-pulse hitting a diffusing screen with 35 ps resolution Devaux95 , thus yielding a 3D imaging. The same technique was later used to image an object embedded in a thick diffusing sample Devaux99 . Althought our setup is actually an improved version of that described in Devaux99 , we point out that our conceptual approach is different from the study of the propagation of a wave front. In fact, in our case the propagation variable is fixed.

Our goal in this article is to demonstrate the potentiality of the optical gating technique to acquire a high resolution space-time map of short, focused WPs in their comoving reference frame. Furthermore, we devise the capability of the technique to reconstruct both the amplitude and the phase of the WP thanks to the coherence of the nonlinear interaction. We propose a method that is based on quadratic type I interaction in a sum-frequency generation scheme either by a non-collinear second harmonic generation or by a collinear sum-frequency scheme. The latter has been used in jose . Here we discuss the first option, showing that if the interaction angle between the two interacting fields is small enough, then a reliable space-time map of the object pulse can be obtained. This can be achieved if the duration of the gate is much smaller than that of the object to be imaged. A holographic interpretation of the method permits to gain insight into the process of up-conversion of the space-time slices of the object into the SF field, and to prove that the coherence of the SF process is able to reconstruct the wavefront in both amplitude and phase. Our results confirm this possibility. The theoretical discussion of the method is followed by section 3, where we present the set-up and the experimentally reconstructed space-time intensity profiles of a parametric spatial soliton excited by a 1 ps light pulse. For our setting, we estimate a mapping resolution of 200 fs in time and about 10$`\mu `$m in space. The features of the technique are presented in section 4, pointing out the limitations that may arise and discussing the possible implementations in each case. In the last section the main conclusions are presented.

## 2 Description of the technique: intensity and field reconstruction

In this section we explicitly show how a non-collinear sum-frequency (SF) scheme can be exploited to get high resolution space-time intensity maps of an unknown light wave packet with a space-time structure (object wave). The recovery of a 3D intensity map is obtained by means of a short reference pulse which provides a time gating inside a nonlinear (NL) crystal, and generates a SF signal containing the information about a set of 2D slices of the object obtained by changing the reference delay. We first discuss the case for the reconstruction of the object intensity profile, and then we show how the intrinsic coherence of the SF process allows in fact for a truly holographic recording of the unknown object.

### 2.1 3D Intensity profile mapping

Let us denote the object ($`\overline{E}_O`$) and reference ($`\overline{E}_R`$) wave packets as follows:

$`\overline{E}_O`$ $`=`$ $`E_O(x,y,z,t)e^{i[\omega _1tk_z(\omega _1)zk_x(\omega _1)x]}+c.c.`$ (1)
$`\overline{E}_R`$ $`=`$ $`E_R(x,y,z,t)e^{i[\omega _2tk_z(\omega _2)z+k_x(\omega _2)x]}+c.c.`$ (2)
where the complex functions $`E_O(x,y,t,z)`$ and $`E_R(x,y,t,z)`$ are the slowly varying envelopes of two waves with frequencies $`\omega _1`$ and $`\omega _2`$. Note that in this form the equations describe two wavepackets propagating in the positive $`z`$ direction and colliding at an angle $`\theta =2\mathrm{arctan}(k_x/k_z)`$ in the $`xz`$ plane (here $`k=\sqrt{k_x^2+k_z^2}=2\pi /\lambda _0`$). For a SF generation process occurring inside a quadratic nonlinear crystal, the polarization source giving rise to the SF can be written as:

$`P_{SF}2E_OE_Re^{i[\omega _3tk_z(\omega _3)z]}+c.c.`$ (3)
where the phase and energy matching conditions $`k_z(\omega _3)=k_z(\omega _1)+k_z(\omega _2)`$, $`k_x(\omega _1)=k_x(\omega _2)`$ and $`\omega _3=\omega _1+\omega _2`$ have been used. The SF field propagates along $`z`$ direction and has a slowly varying envelope function that we will indicate by $`E_{SF}(x,y,t,z)`$. Now we introduce the following assumptions: 1) small depletion of both the $`E_O`$ and $`E_R`$ fields; 2) negligible diffraction and dispersion within the crystal; 3) equal group velocities of the object, reference and sum-frequency pulses, namely $`u_O`$, $`u_R`$ and $`u_{SF}`$, that is $`u_O=u_R=u_{SF}=u`$. Note that all these assumptions approximatively hold as long as the thickness of the crystal $`\mathrm{\Delta }z`$ is small compared to the characteristic lenghts of the system (nonlinear length, dispersion and diffraction lengths, pulse walkoff length). These assumptions allow to find a travelling reference frame for all the propagating pulses by introducing the retarded time $`\tau =tz/u`$. The general partial differential equations describing the interaction process then reduces to an ordinay differential equation for the envelope $`E_{SF}`$. If we also introduce a time delay $`\tau _i`$ on the reference wavefront, the equation takes the form:

$$\frac{dE_{SF}(x,y,\tau ,z)}{dz}=i2\sigma E_O(x,y,\tau ,z)E_R(x,y,\tau \tau _i,z)$$
(4)
where $`\sigma `$ is the nonlinear coupling term. The last equation is readily integrated and, if the mixing crystal is placed at position $`z_0`$, it reads:

$$E_{SF}(x,y,\tau ,z_0)=i\sigma 2E_O(x,y,\tau ,z_0)E_R(x,y,\tau \tau _i,z_0)\mathrm{\Delta }z$$
(5)
A deeper discussion about the meaning of this expression will be given in the following subsection. Here we just point out how the intensity profile of the SF field can be used to retrieve the 3D intensity map of the object. More precisely, for a given lag time $`\tau _i`$, the spatially dependent SF fluence profile (the CCD signal) $`S(x,y,\tau _i,z_0)`$ recorded just at the exit face of the mixing crystal is given by:

$$S(x,y,\tau _i,z_0)(\sigma \mathrm{\Delta }z)^2_{\mathrm{}}^+\mathrm{}|E_O(x,y,\tau ,z_0)|^2|E_R(x,y,\tau \tau _i,z_0)|^2𝑑\tau $$
(6)
This expression provides the convolution between the intensity profiles of the object and reference wave packets, lagged in time by $`\tau _i`$, and shows that the signal recorded by a CCD sensor is a linear function of the intensity of both the object and reference wavepackets. In the particular case in which the reference wave is spatially homogeneous in the transverse $`xy`$ plane, and temporally much shorter than the object, we can write expression 6 in the following form:

$$S(x,y,\tau _i,z_0)(\sigma \mathrm{\Delta }z)^2|E_R|^2I_O(x,y,\tau _i,z_0)$$
(7)
where $`I_O=|E_O|^2`$.

By imposing a set of delays ($`\tau _i`$, $`i=1\mathrm{}n`$) to the reference WP with respect to the object, a reliable 3D reconstruction of the WP structure can be achieved by the collection of the $`n`$ images $`S(x,y,\tau _i,z_0)`$. By changing the plane $`z_0`$, the temporal evolution of the WP can also be obtained.

Notice that the model described above does not take into account the dispersion of the mixing crystal. Therefore, the model as it is predicts no limits for the resolution of the map as long as arbitrary short reference pulses are available. Actually the real mixing crystal has a finite bandwidth that limits the spatiotemporal resolution of the obtainable maps. Therefore a careful evaluation of the dispersion characteristics of the mixing crystal have to be gauged as ultrashort pulses are either investigated or used as a reference. More details are discussed in section 4.

### 2.2 3D field reconstruction

The method described above is not limited to the intensity reconstruction, but can be also implemented for the field reconstruction as can be stated from equation 5. This point can also be explained in terms of a holographic description of the process, thus bringing to a deeper understanding of the imaging reconstruction process.

As stated in Denisyuk99 ; Denisyuk00 ; staselko , for a plane wave reference the generated wavefront through the NL interaction behaves as a conventional hologram recorded at frequency $`\omega `$ and illuminated with a radiation at $`2\omega `$. As a consequence, the position, scale, resolution and all the other properties of transformation of the reconstructed image can be predicted by means of the ordinary laws of holography. This hologram is recorded and reconstructed at a time, and exists only when light propagates inside the crystal. Yet we also point out that we obtain the 3D map of our pulse by collecting a set of independent 2D holograms by slicing the object pulse at different delays. Nevertheless, according to the holographic properties of the SF process, the slicing can be done with the mixing crystal at any $`z`$ from the real object to be recovered and the information recorded is enough to reconstruct the slice.

In order to clarify this point, let’s assume to perform an experiment in which we reconstruct a slice of an object WP (Fig 1.a). At any distance $`z`$ (Fig 1.b) the reconstructed WP will correspond to the propagated version of the one at $`z_0`$. Let us consider the case when the object WP has a bandwidth small enough with respect to the SF bandwidth that our simplified model applies, the pulse does not have appreciable angular dispersion222In the case of angular dispersion, the diffraction drives an effective group velocity dispersion also in the vacuum, then the spatial evolution cannot be separated from the temporal one (see for example sonajalg , zozyula ) and the spatial and the temporal evolution can be separated. Under these assumptions, any slice recorded in the confocal configuration (see Fig. 1.a) could also be reconstructed when the mixing crystal is displaced at position $`z`$ and the field to be converted is the propagated one. This is possible by exploiting the holographic properties, provided that the detecting system is set to reconstruct the virtual image of the slice. This image is placed at a distance 2$`z`$ far from the mixing crystal (see Fig. 1.c), its transverse size being identical to that of the object (see Denisyuk00 , and also note that we are working in the degenerate case, $`\omega _1=\omega _2`$).

The possibility to recover the intensity profile of a virtual image comes from the complete wavefront (field) reconstruction arising from the SF process coherence and contained in Eq. 4)333We suggest for example the possibility to use an interferometric or heterodyne device in order to obtain a complete charaterization of the field. By measuring the field instead of the intensity, a remarkable increment of the dynamical range is also possible..

This important property also allows to get intensity maps with a wider dynamical range. For example, in the cases when the object is so intense so that the undepleted pump approximation does not hold, we can get the WP profile by displacing the NL crystal to a plane where the intensity is reduced and then recover the intensity profile at the desired plane $`z_0`$ by suitably moving the imaging system (Fig 1.c)).

## 3 Experimental results

We prepared several experiments in order to prove the possibility to recover 3D maps of short WPs and to show their holographic properties.The experimental set-up is sketched in Fig. 2.

The 1 ps pulses of a Nd:glass laser source (TWINKLE, Light Conversion, wavelength 1055 nm) are splitted on two lines by means of a beam splitter. In the first line, the laser pulse is frequency doubled in a KDP crystal and focused on a 15mm-long litium triborate (LBO) crystal. For pulses about 1 $`\mu `$J in energy , spatial solitons are formed in this crystal because of the optical parametric generation process DiTrapani98 . The gaussian temporal profile causes the spatial soliton to be formed only in the central part of the pulses where the intensity is higher Barthelemy02 . This leads to a non-trivial space-time structure in the output pump wave packet Minardi03 , which we will consider as the object. Furthermore, the object has been magnified 4 times by means of a two-lens telescope imaging the LBO exit face into the NL mixing crystal ( a 100 $`\mu `$m thick, $`\beta `$ barium borate crystal, BBO). The beam expansion has been necessary to: i) reduce the beam intensity and therefore fulfill the undepleted mixing requirement; ii) to avoid information loss related to the finite spatial resolution of the mixing process. In the second line, a 200-fs reference pulse is produced at the wavelength of 527.5 nm by means of a second-harmonic pulse compressor Stabinis91  and expanded 3 times by means of a telescope. Both lines are recombined in the thin BBO crystal, cut and oriented to generate the non-collinear SF from the object and reference beams. An external incidence angle of $`6.5^{}`$ between the propagation directions has been chosen. The delay between the two pulses can be varied by means of a suitable delay line placed on the reference pulse line. The non-collinear SF radiation is spatially selected by an aperture, then selected in frequency by means of coloured filters and the plane of the BBO is imaged onto the CCD sensor (PULNIX TM6CN). Finally, both the BBO crystal and the imaging system could move independently on a rail along the optical axis.

At first we focused the object inside the BBO crystal and adjusted the position of the imaging system as in Fig. 1.a). By recording the spatial profile of the SF radiation as a function of the reference pulse time lag (steps of $`\mathrm{\Delta }\tau _i=66.6`$ fs has been used, where $`\mathrm{\Delta }\tau _i=\tau _i\tau _{i1}`$), we have retrieved the 3D isointensity maps of the object pulse at the exit face of the LBO crystal.

Fig. 3.a shows three different intensity levels of the pulse in the $`(x,y,\tau )`$ space, while the corresponding contour plot of the $`(x,\tau )`$ plane section is depicted in Fig. 3.b. The space-time maps clearly show that the pulse is formed by a spatially focused structure followed by a diffracted tail. The accuracy of this reconstruction has been successfully tested by comparing the experimental plots with the results of a 3D-numerical simulation of the formation process of the object pulse (see Minardi03 ).

In order to check for the validity of the holographic interpretation in this process, we have reconstructed the same WP of Figure 3 by moving the NL crystal a distance z far from the previous position and by scanning the position of the imaging system (CCD+lens) along the SF propagation direction in order to find out the position of the (virtual) reconstructed image. First we fixed the time lag between the object and the reference and selected a slice of the object corresponding to a narrow focused spot of $`40\mu `$m in diameter inside the BBO crystal (corresponding to case a) in Fig. 1). Furthermore, the displacements $`z`$ of the BBO crystal have been chosen large enough to ensure that the propagated wavefronts had lost the transverse structure. As a rough estimate, for the focused part of the object, the Rayleigh range is about 7 mm long, while we spanned distances from -15 mm to 15 mm. In Fig. 4, data show the position of the virtual image against the position of the real object from the crystal (see Fig. 1.c)), showing a remarkable fit to a straight line which slope is close to the value of 2, according to what discussed in the precedent section. The error bars shown in the figure indicate the estimated uncertainties in the reconstructed image plane position, measured by scanning the whole Rayleigh range and by extracting the position of the central point. In the same figure the two insets show two isointensity maps recovered for the case when the imaging system is focusing directly the BBO plane (z=0) and when the BBO crystal has been moved a distance z=15 mm (the maximum propagation distance we imposed to the BBO). The agreement between the two intensity profiles proves the reliability of our method to work with virtual images in recovering 2D slices of WPs like those we used here.

In order to verify that we have always been operating in the undepleted regime, we have measured the dependence of the generated SF field as a function of the object-pulse energy at $`z=0`$. The results are presented in Fig. 5, where the peak fluence (as registered from the CCD images) is plotted vs the transmission of the neutral-density filters that attenuate the object. The costant slope confirms the absence of any saturation in the SF process. Data also indicate that a slight overestimate of the background has been done during the profile acquisition (zero SF field is found for 15 % filter transmission).

## 4 Features of the technique

In the previous sections we have proved that the described technique is able to perform the 3D mapping of objects similar to those we used for our measurements. Anyway, one can expect that in general the reliability of the technique to recover the slices could be affected by particular effects dealing with the interaction geometry, the mixing crystal, and the features of the interacting pulses, in particular for broad-band, ultra-short and chirped pulses. Therefore we discuss here three main features of our technique, namely, the duplication bandwidth of the mixing crystal, the interaction angle and the reference pulse shape, that bring to devise some possible limitations to the fidelity of the technique.

### 4.1 Space-time resolution

As we briefly mentioned above, the interaction model we describe does not take into account the actual finite bandwidth of the mixing crystal, limiting the spatiotemporal resolution of the maps. In fact, in our experiment the BBO crystal has been chosen thin enough that the converted bandwidth was large compared to the object one.

As a matter of fact, the spatio–temporal resolution of the slices is dictated by the maximum range of angles and frequencies over which the conversion is effective (for very broad–band objects the angular blurring that could arise from the phase matching condition should be considered). The key parameter here is the maximum phase mismatch between the fundamental and the SF waves at which the conversion efficiency vanishes, $`\mathrm{\Delta }k_{max}`$. On the basis of the existing literature niko , we have estimated the maximum bandwidth converted by the crystal as the FWHM of the efficiency curve, yielding to a temporal bandwidth of about $`630`$ cm-1 and an angular bandwidth of about $`510`$ cm-1. By considering only the angles for which the conversion efficiency is higher than 1/e of the maximum value devaux , we obtain the remarkable resolution of details approximately 10 $`\mu `$m in size (100 lines/mm). We can notice that the holographic interpretation of the wavefront reconstruction leads to an easy description of the imaging process. Furthermore, our optical system could be used in principle to achieve a microscopy of ultrashort objects shorter than the one we used here jose , although the limit of very shot WP is detrimental for the holographic reconstruction of the phases.

### 4.2 Interaction angle

A point to be discussed concerns the influence of the noncollinear geometry and pulse chirp on the reconstruction of the hologram. In principle the holographic interpretation is strictly valid for collinear geometry only (when $`k_2=2k_1`$) Denisyuk00 . Although this condition is not strictly fulfilled in our experiment, we have maintained the interaction angle small enough to make this disturbance negligible. However, a temporal chirp in the object pulse could give rise to a spatial phase distortion of the SF wavefront Stabinis91b . Although this does not affect the results obtained with the confocal configuration (see Fig. 1.a), it could distort the maps obtained in the holographic one (Fig. 1.c). We expect that this effect is really rampant only when complicated ultrashort pulses with strong chirp are considered, or large interaction angle are employed. We are confident that our system was operated far from this limit, since, as we checked from the data, our experimental holograhic maps do not show any appreciable enhancement of the astigmatism.

As it is well known, the non-collinear SF scheme is largely used to get single-shot autocorrelation traces of ultrashort pulses Gyuzalian79 . In fact, because of the geometry of the interaction, it is easy to show that the time of the interaction between the object and reference pulses depends on the transverse coordinate of the intersecting planes (see (zewail, , pag. 426-428)). Therefore, if the interacting angle is large, we expect that distortions of the spatial profile along this direction could affect the space-time maps. However, it turns out that in the confocal configuration (i.e. as in Fig. 1.a) the effect of the interaction angle is merely that of skewing the intensity map, so that the actual time-axis of the object pulse is not parallel to that measured on the delay line. This is evident in the map depicted in Fig. 6, obtained with an external crossing angle of $`16^{}`$.

In the paraxial approximation, the slope of the object-pulse time-axis on the space-time map is given by:

$$\mathrm{tan}\beta =\frac{c}{n}\mathrm{sin}\frac{\gamma }{2}$$
(8)
where $`\gamma `$ is the angle between the propagation directions of the object and the reference wavepackets inside the crystal, and $`n`$ is the refractive index of the medium. It is easy to recognize this formula as the calibration relation of the non-collinear, single-shot intensity autocorrelator Gyuzalian79 ; zewail .

### 4.3 Reference pulse shape

Finally the last point to be discussed concerns the shape of the reference pulse. In our experiments it was only slightly affected by the presence of satellites because of the generation via the pulse compressor Stabinis91 , the amplitude of these structures being small enough to be considered negligible to our aim. Actually, the recovered intensity profile of the object would be distorted in case the reference pulse has a structure more complicated than the single peak envelope, as Eq. 6 clearly points out. This suggests that a careful measurement of the reference shape (for example carried out by means of an autocorrelation technique) can be used to retrieve the real intensity profile by means of a deconvolution procedure. We point out, however, that any deconvolution unavoidably introduces an extra spurious noise which could degrade the final quality of the mapping deconvolution ).

## 5 Conclusions

We have shown that the 3D intensity maps of optical WPs with a complex space-time structure can be retrived by an optical gating technique. The method allows the reconstruction of the WP in its comoving reference frame and, by exploiting the holographic properties of a slightly non-collinear degenerate, sum-frequency process, we have shown that a complete amplitude and phase reconstruction is actually obtained. The maps are obtained by suitably imaging the second harmonic radiation obtained by cross-correlating an object pulse with a much shorter plane wave-packet delayed in time. The holographic properties of the generated radiation have been carefully tested experimentally, and the distorsions of the intensity maps introduced by the non-collinear geometry have been discussed in detail. Finally, theoretical considerations point out that the ultimate resolution of our set-up is in the order of 100 lines/mm in space, and of a few fs in time. However, the choice of the reference pulse limits the actual time resolution to about 200 fs.

We forsee that the developed technique will benefit all the fields where a space-time mapping of light pulses is relevant, such as the investigations on the reshaping of ultrashort pulses propagating in non-linear materials Lantz02a ; Lantz02b ; PDT00 ; Minardi03 ; jose . Moreover, the holographic features of the technique might be exploited to fully reconstruct the field of an object pulse and its evolution during the propagation.

This work was partially supported by MIUR (COFIN01 and FIRB01), the European Commission EC-CEBIOLA project (ICA1-CT-2000-70027) and DGI BFM2002-04369-C04-03 (Spain). The work of J.T. is supported by a postdoctoral grant from Ministerio de Educacion Cultura y Deporte (Spain).

$``$ Permanent address, Department of Physics and Nuclear Engineering, Universitat Politécnica Catalunya, 08222 Terrassa, Spain.