1. INTRODUCTION

Nonequilibrium pattern formation occurs in many physical systems, giving rise to the emergence of order on macroscopic scales [1,2]. The theory of hydrodynamics is a paradigm for understanding these phenomena, and is applicable to many branches of physics, chemistry, biology, astrophysics, and other sciences [3,4]. This is often applied to many body systems subject to nonlinear coupling in a dissipative environment with external fluxes. In physics, one of the most studied hydrodynamic effects is the theory of fluid flows near the Rayleigh–Bénard instability [1,2]. The simplest nontrivial model for fluid convection that displays pattern formation due to convective instabilities was derived by Swift and Hohenberg [5], where nonlinear coupling of fluctuations was included to demonstrate the failure of mean field theory for critical exponents.

Here, we treat pattern formation and universality in a paradigmatic nonequilibrium quantum system: the nondegenerate parametric oscillator (OPO). Theories of a similar nature have been applied to nonequilibrium spatially extended structures in lasers and other related systems, with an emphasis on universal behavior of phase transitions, pattern formation, and self-organization [6,7]. Yet downconversion in a nondegenerate parametric system can display new possibilities not found in these simpler cases. In particular, one can have entanglement and Einstein–Podolsky–Rosen (EPR) paradoxes [8–11].

In the present paper, we investigate a new type of critical point phase transition in this nonequilibrium quantum system in order to understand the universality class. In a type II OPO, there are two downconverted fields with orthogonal polarization. Hence, the order parameter is a complex or vector field in two dimensions. The two components of this vector field are associated with the polarization degrees of freedom of the downconverted radiation field. This is a quantum system driven to a phase transition far from thermal equilibrium. It is also known to display strong quantum entanglement and EPR correlations [8–10] in the case where there are two correlated output modes. We wish to understand the behavior of this phase transition using a first-principles analysis [12,13] of the relevant master equation.

Experimentally, the OPO is now a mature technology with both commercial and fundamental applications. Following initial theoretical predictions [14–16], the quantum limited type I OPO was investigated experimentally by Wu et al. [17], demonstrating quantum squeezing. Later the type II case, in a triply resonant cavity, was used to experimentally demonstrate continuous variable EPR correlations [11], also originally predicted theoretically [12,13]. These initial experimental investigations were in few-mode devices. Both below- and above-threshold experiments have been carried out, as close as $\pm 1\%$ of the critical point [18–23], confirming predictions of thresholds and conversion efficiency. Operation at the critical point results in low-frequency critical fluctuations and non-Gaussian behavior [24,25].

Spatially extended pattern formation in a degenerate or type I OPO has also been analyzed previously [26]. It is related to the Lifshitz phase transition [27]. This is a model used to describe the phase transition to a modulated magnetic phase [28–30]. In this simplest case one has a two-dimensional, planar system with a scalar-order parameter [31], which has known universality properties. A Swift–Hohenberg equation was derived for spatially extended nondegenerate type I OPOs with flat end mirrors [32], but ignoring fluctuations. OPO experiments have been recently extended to these types of multimode devices [33–35], with type II as well as type I parametric downconversion. The experimental situation is that while multimode experiments have mostly used confocal mirrors, the first type II multimode planar mirror experiments have been carried out [35].

Studies of the Swift–Hohenberg equation in the neighborhood of the critical point, the Lifshitz point, as well as pattern formation for lasers have been discussed [36]. The Lifshitz point is similar to the tricritical point [37]. Tricritical points occur in different physical systems. They correspond to a point in a phase diagram where two lines of ordinary critical points meet and terminate [38–40]. For critical points and tricritical points there are two critical dimensions. The upper critical dimension refers to the one above which the critical exponents have classical values [37]. The lower critical dimension refers to the smallest dimension for which there is a true phase transition, due to increasing fluctuations as the dimension is decreased.

Systems like these can entangle large numbers of modes, with quantum entanglement increasing as threshold is approached. Currently, experiments near the critical point are sensitive to classical fluctuations [24,25] and heating effects [35]. With improved stabilization methods, we expect that these technical problems can be overcome. However, noise due to quantum fluctuation effects will remain.

Theoretically, the usual approach for two-dimensional nonequilibrium problems is the Landau–Ginzburg (LG) equation [41,42], which was first used for understanding superconductivity near threshold. Normally, the LG equation is derived using approximations like adiabatic elimination, and is taken as an effective equation for the physical system. There is now an increased interest in extended spatial and multimode structures [33,34,43–46] in the quantum optical parametric oscillator (OPO) [17,47,48], due to availability of multiple transverse mode cavities [49] and quantum imaging control [50]. It is important in these cases to understand how quantum noise enters the dynamical equations [51].

We show that planar type II downconversion creates a nonequilibrium system with a similar general type of symmetry to the Berezinskii [52] and Kosterlitz–Thouless (BKT) models [53], but with an isotropic Lifshitz point, first studied in magnetic systems by Hornreich et al. [37]. Physically, this is not a classical fluid or magnetic system, but rather is a quantum system, driven into a nonequilibrium critical point [54,55] far from thermal equilibrium. The nondegenerate planar OPO is fundamentally different to the usual BKT model, having a quartic rather than quadratic momentum dependence in the linear response function. This places it in a similar category to the Swift–Hohenberg model and next-nearest neighbor lattice models. Our model can also display strong EPR entanglement and other nonclassical properties in addition to the Lifshitz point behavior. The parametric model therefore provides a novel path to the investigation of unusual classical and quantum noise effects.

Phase transitions with this general type of symmetry are continuous yet break no symmetries. There is also no ferromagnetism in this system, as this is prohibited by the Mermin–Wagner theorem [56]. Berezinskii predicted a new type of phase with correlations that decay slowly with distance with a power law. The phase transition for a ferromagnetic phase is prevented by the appearance of vortex and antivortex pairs. Many quantum phase transitions in two dimensions belong to this class. The continuous version of the XY or Ising models are often used to model systems that possess order parameters with a symmetry of this type, e.g., superfluid helium, liquid crystals, two-dimensional Bose–Einstein condensate (BEC), and others. While having the same order parameter and space dimension as this case, we show that the type II OPO has a different universality class.

Our approach is to start from the usual master equation model that describes a type II OPO with transverse modes. Here we consider that the OPO is nondegenerate in polarization. We then map the coupled Heisenberg equations into the positive P-representation [57] of the density matrix. This choice is made because the positive P-representation allows us to exactly map the density matrix evolution into an equivalent stochastic equation.

We organize this paper in the following way. First we describe the Hamiltonian model and derive the equation of motion for this open system. Next, we use the positive P-representation for mapping these equations into a Langevin type, which can either be treated numerically or via analytic approximations. We will mostly describe the unsqueezed quadratures of the field, which have a large similarity with the corresponding magnetic system. At this point we can recognize the universality class. We give both a Gaussian approximation to the correlation functions of the unsqueezed quadratures, and high-precision numerical simulations of the non-Gaussian corrections. We show that, even though the system is far below its upper critical dimension, the non-Gaussian character of the fluctuations is relatively small, and the intensity fluctuations are nearly factorizable.

In a following paper we will focus our attention on the squeezed field quadratures, to give a spatial map of quantum squeezing and EPR entanglement.

2. MODEL

The system of interest comprises an optically driven planar Fabry–Perot cavity or interferometer with a nonlinear medium that possess a parametric nonlinearity. The nonlinear crystal is cut to give a type II phase matching that couples a pump field to two downconverted fields having an orthogonal polarization. The cavity is pumped with a spatially extended coherent light with frequency ${\omega }_0$, with a transverse spatial profile. In the simplest case we consider this pump to be a plane wave.

A. Quantum Hamiltonian

The outgoing downconverted light, amplified inside the cavity, develops structures and patterns due to diffraction, nonlinear coupling, and detuning between the wavelength of the downconverted field and the cavity size. This depends on the modal decomposition of this cavity. The mirrors have parameters that can be controlled by experimentalists. The tunable parameters include the reflection coefficient for each mode, and the cavity detunings for each mode.

Our model for this system is similar to many earlier treatments of driven nonlinear optical cavities [12,14,15,43–46,58–63]. It includes a linear coupling between the external electromagnetic field modes and the internal cavity modes, owing to a partially transmitting mirror. The cavity and mirror parameters determine both the coupling to the driving field and the decay rate of the cavity or interferometer. We note that for a low-$Q$ device, it is important to use nonorthogonal quasi-modes [59]. Here, we assume the opposite case of a thin, high-$Q$ extended planar cavity, so that the external modes are simply plane-wave modes.

The quantum Hamiltonian in the interaction picture has four main terms that can be summarized by the following expression,

The free-evolution Hamiltonian that accounts for diffraction inside the planar cavity is

The interaction Hamiltonian representing the coupled modes inside a crystal with ${\chi }^{(2)}$ nonlinearity is given by [65]

From now on, for definiteness, labels 1 and 2 stand for polarizations, and we shall focus on the degenerate frequency case with nondegenerate polarization. For dimensional reasons, $\chi \propto {\chi }^{(2)}/\sqrt{\ell }$, where ${\chi }^{(2)}$ is the Bloembergen nonlinear polarizability coefficient and $\ell$ is the intracavity longitudinal mirror spacing.

Following standard input–output theory derivations [14–16,60–63], the Hamiltonian term associated with the input laser pumping in a rotating frame at frequency ${\omega }_L$ is

B. Master Equation

The nonunitary evolution of the system comes from the coupling between the cavity modes and the output modes. This can be treated as a quantum Markovian process that simulates a bath interaction. We carry out this calculation in a type of interaction picture so that the interaction picture operators evolve according to a reference Hamiltonian ${\widehat{H}}_0$, given by

We note that, although we focus on the master equation approach here, it is sometimes useful to write the time evolution in a complementary quantum Langevin formulation. This formulation is given in Appendix A.

3. STOCHASTIC EQUATIONS IN THE POSITIVE P-REPRESENTATION

As nonlinear operator equations are not generally soluble, it is more manageable to map the operator equations into c-number form. In this approach, a master equation is transformed into a positive-definite Fokker–Planck equation using operator identities that map the operator terms in the master equation into differential operators [60,63,66]. To do this we have to use a phase-space representation of the master equation.

Phase-space representations in a classical phase space do not give a positive-definite equation. An example is the Wigner representation [67]. While this is exact, it is not able to be mapped into Langevin equations, unless one either truncates or uses higher-order noise [68]. One can also linearize the Hamiltonian and obtain an approximate Wigner diffusion [54,55]. Another approach is the Husimi $Q$-function [69], where, in order to obtain a positive Fokker–Planck equation, an unphysical constraint on the phase-space trajectories has to be used [70].

Here we wish to have the ability to treat nonequilibrium structures without restrictions. For this purpose, the most useful representation is the positive P-representation [57], which is an extension of the Glauber–Sudarshan P-representation [71,72] into a phase space of double the classical dimensions. Unlike the P-representation, which is singular for nonclassical states, the positive P-representation is well defined, positive, and nonsingular for any quantum state. This approach allows us to map the density matrix equation into a Fokker–Planck equation on a nonclassical phase space. In the positive P-representation stochastic averages give normally ordered quantum expectation values. A brief description is given in Appendix B.

The stochastic field partial differential equations are given by an extension of our earlier work [31]:

The stochastic fields ${\xi }_k$ that describe the quantum noise are complex and Gaussian, whose nonvanishing correlations are

A. Critical Driving Field

As a first investigation, we treat the classical approximation, which has also been analyzed in some earlier work [32,43–46,73–81]. Here one assumes that all noise is negligible, so that $A_i^{+}=A_i^{*}$, which gives equations in the form

Solutions with negative intensities are unphysical. There is a positive, above-threshold solution if $|\chi \mathcal{E}|>|z|$, which gives an identical critical field to Eq. (15). For $|\mathcal{E}|>{\mathcal{E}}_c$ there is a transfer of energy from the pump to the signal and idler modes, which develop a finite mean intensity. It is the vicinity and just above this critical point that is the main regime of interest in this paper. We note that above threshold, further instabilities exist, including limit cycles and spatial pattern formation [75].

4. ADIABATIC ELIMINATION OF THE PUMP MODE

We now return to the full quantum behavior given by the stochastic equations obtained above. One limit that has an especially simple behavior is found in the case of a rapidly decaying pump mode. We can treat this by means of an adiabatic elimination procedure. Assuming that ${\tilde{\gamma }}_0\gg {\tilde{\gamma }}_1\simeq {\tilde{\gamma }}_2$, and that $\mathcal{E}$ is spatially uniform (that is, we are neglecting pump diffraction), we can perform an adiabatic elimination by using the stationary solution for the pump mode, so that

A. Signal and Idler Equations

The resulting equations for the downconverted modes—often called the signal and idler equations—are, for $i=1,2$,

Although these assumptions simplify the algebra, they are not essential for our main conclusions, which mostly rely simply on the fact that we now have a two-dimensional order parameter rather than a one-dimensional order parameter as found in the degenerate case. We should note that a nonzero pump detuning can excite another mode different from ${\omega }_0$. This could give rise to nonlinear phenomena like bistabilities, nonlinear resonances, and subcritical bifurcation, especially for the case of large detuning [32]. In the case where the diffraction terms are different, there will be a new term which is proportional to the difference of the two diffraction terms, as has been studied elsewhere [32]. This term will be present in the equations even for the case of zero detuning. Since we are interested in the universal behavior of the system we do not include these cases, but we point out that they can give rise to nonlinear behavior.

B. Dimensionless Form

Equation (20) will now be transformed into a dimensionless form, which allows comparisons with other types of phase transitions. First, we define the dimensionless variables $\tau =t/t_0$, $\mathit{r}=\mathit{x}/x_0$, with a scaled Laplacian ${\nabla }_{\mathbf{r}}^2={\partial }^2/\partial r_1^2+{\partial }^2/\partial r_2^2$. It is useful to also define a dimensionless field ${\alpha }_i=x_0{A}_i$ as well. This has an intuitive interpretation as the coherent amplitude in real space, defined relative to a physical area of $x_0^2$. After this transformation, one obtains

5. STABILITY PROPERTIES AND QUADRATURE EQUATIONS

We now wish to transform these equations into quadrature equations that are simpler to investigate. There are very different stability properties for the orthogonal quadratures near the critical point. To investigate this, as a first approximation, we will ignore noise and nonlinear terms of order $\sqrt{g}$ and smaller. The stability of the equations near ${\alpha }_i=0$, to leading order in $g$, is

A. Quadrature Field Variables

To understand the behavior of Eq. (28) in the neighborhood of the critical point, we define complex, dimensionless scaled quadrature fields that are proportional to the critical eigenvectors, as follows [31,54,55]:

Next, we consider the case $\mu \approx 1\gg g^2|{\alpha }_1{\alpha }_2|$, to simplify the noise term. As this is both relatively small and nearly constant in the neighborhood of the critical point, the resulting terms are of higher order in $g$ than the leading terms we wish to include. The resulting equations for these quadratures in the positive P-representation near the critical point are

At this stage, we can make the following remarks. The stochastic quadrature fields $X,X^{+}$ are both complex fields, so they have four degrees of freedom between them, and similarly for $Y,Y^{+}$. They have a correspondence with non-Hermitian operator fields $\widehat{X},{\widehat{X}}^{†}$ and $\widehat{Y},{\widehat{Y}}^{†}$. In general, $Y,Y^{+}$ are not complex conjugate except in the mean, and neither are $X,X^{+}$, since they are driven by the $Y,Y^{+}$ fields. However, as we will show, this picture simplifies when one considers an expansion near the critical point.

B. Critical Point Adiabatic Elimination

We can now perform a second type of adiabatic elimination, which is valid in the neighborhood of the critical point. This takes into account the fact that the fluctuations in the $X$ quadrature become very slow near threshold, while the $Y$ quadrature still responds on the fast relative time scale $1/\gamma$. Formally, we can drop terms of $\mathcal{O}(\sqrt{g})$ where $g\ll 1$, and approximate Eq. (32) as follows [31],

C. Vector Swift–Hohenberg Equations

We will now show that these near-threshold equations are actually coupled or vector Swift–Hohenberg equations [3,4] that represent the leading-order dynamics near threshold of the downconverted modes with the same frequency but orthogonal polarization. To demonstrate this, it is convenient to make the following change of variables:

 

Fig. 1. Dimensionless intensity versus dimensionless pump $\mu$ in the vicinity of the critical point. The point $\mu =1$ corresponds to ${\eta }_1={\eta }_2=0$ and ${\eta }_3=\frac{1}{2}$. Here we have used the parameter $g=0.01$. Results were obtained from a simulation with 300 samples on a $50\times 20\times 20$ numerical grid of $30000\times 96\times 96$ points.

Download Full Size | PDF

Direct simulations were employed, using both a central partial difference algorithm in the interaction picture [83,84] and a fourth-order interaction picture Runge–Kutta method. Two public domain software packages were used and checked against each other [85,86]. This figure is obtained by solving the stochastic equations given in Eq. (43) at $\mu =0.9$, equilibrating for $t=50$ units then scanning the driving field $\mu$ adiabatically (using $\mu =0.9+0.004t$) through the critical point until $\mu =1.1$. There is a rate-independent critical region at the threshold value of $\mu =1$, where the transition is smooth rather than discontinuous, as it would be classically.

In the case of a one-dimensional, real order parameter, this equation was first derived by Swift and Hohenberg [3–5,87,88]. Here it appears as a two-dimensional vector equation, since the order parameter is two-dimensional. This was used to explain the convective roll patterns generated by the Rayleigh–Bénard instability, where the order parameter is the vertical fluid velocity. The real order parameter case is also similar to the Ising model for magnets in two dimensions with next-nearest neighbor interactions [89], where competition between the nearest and next-nearest interactions generates a magnetic modulated phase called the Lifshitz phase [29,30].

Higher-dimensional or complex order parameters, as in the present analysis, are described using a generalized Landau-Ginzburg-Wilson Hamiltonian [37], where the authors also introduce the Lifshitz point. For the present case, the upper critical dimension for classical behavior is at $d=5$, and the classical location of the Lifshitz point is at ${\eta }_1={\eta }_2=0$. Low-dimensional cases should have enhanced fluctuations, without spontaneous magnetization or symmetry breaking, as expected from the Mermin–Wagner theorem [56]. This applies to Heisenberg-type models with higher-dimensional order parameters in two dimensions. As it is valid for general, finite-range interactions, it also holds in our case.

6. CORRELATIONS

All physical quantities that we wish to understand in the double adiabatic limit treated in the previous section come from the solution of Eq. (43). This allows us to calculate the expectation value of any other observable in the vicinity of the critical point. Another way to obtain expectation values is to write the functional probability as a solution of the master equation. Below we develop both methods.

A. Stationary Solution of the Fokker–Planck Equation

We start with Eq. (43) and note that it is possible to write a functional Fokker–Planck equation for the probability density $P(\mathit{X},\tau )$,

B. Stochastic Moments in the Gaussian Approximation

In order to evaluate the moments and spatial correlations we will approximate the nonlinear terms of Eq. (43) using a Gaussian approximation together with a Green’s function approach. Hence we can write Eq. (43) as follows:

C. Lifshitz Point

We now consider the line of points where ${\eta }_2=0$ and ${\eta }_3=\frac{1}{2}$. These are the points we have defined as corresponding physically to zero detuning, with a pump in the vicinity of the $\mu =1$. In this case, the solution for $\tilde{\mathit{X}}(\mathrm{k⃗},t)$ is given by

 

Fig. 2. Dimensionless intensity versus detuning $\mathrm{\Delta }$. The Lifshitz point corresponds to zero detuning. Here we have used $\mu =1$ and $g=0.01$. Results were obtained from a simulation with 60 samples on a $200\times 50\times 50$ numerical grid of $15000\times 96\times 96$ points.

Download Full Size | PDF

In order to obtain this figure, we solve the stochastic equations given in Eq. (43), and scan the detuning which is proportional to the parameter ${\eta }_2$ defined in Eq. (39), so that $\mathrm{\Delta }=-0.5+0.005t$. We notice that the fluctuations depend on the detuning. For a positive detuning there is a decrease of the fluctuations while for a negative detuning the fluctuations increase. A classical Swift–Hohenberg equation with a complex order parameter and nonzero detuning has been treated [92,93], as has the hyperbolic complex Swift–Hohenberg equation [94] and other related studies [95,96]. For negative detuning there is a ring with strong fluctuations, shown in Figs. 3 and 4.

 

Fig. 3. Dimensionless intensity versus detuning $\mathrm{\Delta }$ and transverse momentum $k_x$. For negative detuning there are fluctuations for lower momentum values. Results were obtained from a simulation with 800 samples on a $100\times 50\times 50$ numerical grid of $15000\times 96\times 96$ points. Here we have used $\mu =1$ and $g=0.01$.

Download Full Size | PDF

 

Fig. 4. Left: Dimensionless intensity versus transverse momentum $\mathbf{k}$ for a fixed value of the detuning $\mathrm{\Delta }=0.5$. Here we have used $\mu =1$ and $g=0.01$. A ring is formed for lower values of the momentum. The fluctuations (peaks) are not symmetric due to noise. Right: Dimensionless intensity versus transverse momentum $k_x$ for a fixed value of the detuning $\mathrm{\Delta }=0.45$, and $k_y=0$. Other parameters are as in Fig. 3.

Download Full Size | PDF

D. Spatial Correlations in the Gaussian Approximation

From Eq. (53), the Gaussian correlation function in the momentum space (far field), in the stationary regime is therefore given by

E. Non-Gaussian Behavior and Universality

In order to verify these analytic results and investigate non-Gaussian correlations, numerical simulations were carried out of the original stochastic partial differential equations of Eq. (43). Small time steps are needed to treat the quartic growth of the squared Laplacian term in momentum, together with large sample numbers to obtain a low sampling error. In initial investigations, we used a $10\times 10\times 10$ numerical grid of $8000\times 64\times 64$ points in $t$, $x$, $y$, respectively. Employing a fine numerical grid of $16000\times 64\times 64$ points to check convergence in time step, the final steady-state correlation result converged to $〈\mathit{X}·\mathit{X}〉=0.264\pm 0.005$. This was close to the Gaussian value, but with a relatively large sampling error.

In order to understand the quantitative difference between the exact and Gaussian results, a more precise differencing technique was used. This variance reduction or differencing technique simulates the difference between the full sample path and the Gaussian approximation [54,55]. The results were in agreement with a direct simulation, but gave much more rapid convergence. This was carried out as follows. First a mean field variable, $\tilde{\mathit{X}}$ was simulated, in the Gaussian approximation, using

 

Fig. 5. Growth of non-Gaussian correlations, $〈{|X|}^2-{|\tilde{X}|}^2〉$ versus time $\tau$, starting from $\mathit{X}=0$. Error bars were obtained from comparing a coarse (5000 step) and fine (10000 step) simulation. The sampling error is indicated by the upper and lower solid lines, with a standard deviation of $\pm 0.00025$. This gives the steady-state correlation result $〈\mathit{X}·\mathit{X}〉=〈{|X|}^2〉=0.2574\pm 0.0003$. Results were obtained from a simulation with 3200 samples on a $10\times 20\times 20$ fine numerical grid of $10000\times 48\times 48$ points.

Download Full Size | PDF

In summary, the full statistical calculation gives increased critical fluctuations due to non-Gaussian effects, but this increase is relatively small. At large transverse momentum there is no measurable difference between the Gaussian and exact results, shown in Fig. 6. The deviation from the Gaussian approximation vanishes rapidly as higher-order transverse momenta are investigated. Small values of this difference are observed only at zero transverse direction.

 

Fig. 6. Left: Steady-state momentum correlations, $〈{|X(\mathbf{k})|}^2〉$ versus transverse momentum $\mathbf{k}$, starting from $X=0$. Right: Steady-state non-Gaussian correlations in momentum space, $\mathrm{\Delta }\ln 〈{|X(\mathbf{k})|}^2〉=\ln 〈{|X(\mathbf{k})|}^2〉-\ln 〈{|\tilde{X}(\mathbf{k})|}^2〉$ versus momentum $\mathbf{k}$. Here we consider ${\eta }_1={\eta }_2=0$ and ${\eta }_3=\frac{1}{2}$. Other parameters are as in Fig. 5.

Download Full Size | PDF

7. CONCLUSION

We have shown that parametric downconversion in a type II parametric planar cavity leads to a Swift and Hohenberg type of stochastic equation for the leading terms in the critical fluctuations, but with a vector order parameter. This combines the rotationally invariant symmetry properties of the X–Y model with the higher-order Laplacian of a Lifshitz magnetic phase transition. Surprisingly, these fluctuations are not thermal in origin, but come instead from the quantum fluctuations associated with parametric amplification.

This model can be approximately treated for the critical fluctuations with a Gaussian factorization. However, a careful numerical treatment shows that non-Gaussian critical fluctuations occur. These are responsible for enhanced intensity correlations, but are reduced at large transverse momenta due to the momentum dependence of the linear propagator. As techniques improve, we expect that this novel, nonequilibrium critical point will become accessible to experimental studies.