Phase transitions in an Ising chain interacting with a single mode cavity field

We investigate the thermodynamics of a combined Dicke- and Ising-model which exhibits a rich phenomenology arising from the second order and quantum phase transitions from the respective models. The partition function is calculated using mean field theory, and the free energy is analyzed in detail to determine the complete phase diagram for the system. The analysis reveals both first- and second-order Dicke phase transitions into a super-radiant state, and the cavity mean-field in this regime acts as an effective magnetic field, which restricts the Ising chain dynamics to parameter ranges away from the Ising phase transition. Physical systems with a first order phase transitions are natural candidates for metrology and calibration purposes, and we apply filter theory to show that the sensitivity of the physical system to temperature and external fields reaches the 1/N Heisenberg limit.


Introduction
The understanding of the remarkable and useful properties of matter in different phases and of the critical behaviour near phase transitions presents ongoing challenges in theoretical and experimental physics [1,2,3,4,5]. Entirely new types of phase transition phenomena become relevant with the ability to control and engineer microscopic interactions and systems, e.g., in cold atom experiments with spinor Bose-Einstein condensates and with Fermi gases [6,7]. Indeed, a whole branch of quantum computing research attempts to use candidate systems for quantum computing as quantum simulators in which quantum gates are operated so as to simulate suitable inter-particle interactions [8] and in this way implement theoretical phase transition models in a quantum analog computer.
Local interactions between nearest neighbours in a spin chain lead to Ising-, Heisenberg-, and other interaction models with phase transitions at definite values of the interaction strengths and external controllable parameters, such as a bias magnetic field. While these interactions are reliable models of, e.g., magnetic interactions in solids, they can also be engineered exactly among trapped atoms or ions, with the added experimental possibility to control the sign and magnitude of the interactions with laser beams and the spin temperature by optical pumping [9].
Atomic and optical systems also permit the engineering of interactions between a large number of atoms and a quantized oscillator mode. Such systems are implemented in various schemes for quantum computing, where the oscillator mode is used as a data bus between the atomic quantum bits. Beyond a critical coupling strength to the oscillator and below a critical temperature the system undergoes a phase transition, and the thermodynamic ground state acquires a macroscopic excitation of the oscillator mode. This phase-transition was first discussed by Dicke [10], and the Dicke phase transition has since then been studied extensively [11,12], and was recently observed in experiments with cold atoms in an optical cavity [13] using techniques similar to those described in [14].
It is conceivable that one can implement both the Ising and the Dicke interaction Hamiltonians in many different ways using atoms and cavities, nano-mechanical devices or collective vibrational motion of the atoms or ions. The partition function of such a combined Dicke-Ising model was determined recently in the thermodynamic limit [15], and in this article we study the phase diagram and further properties of the system and possible applications.
First order phase transitions represent discontinuous changes and hence a high sensitivity to variations in the values of the physical parameters of the model near the critical point. Sensitivity of quantum systems in fundamental metrology is a very active research field where much research has been devoted to identify how the sensitivity depends on the number of particles N . The "standard limit" 1/ √ N , is replaced by the "Heisenberg limit" 1/N within different realizations of non-interacting particles [16,17,18,19], while more rapid decrease with particle number has been proposed in different models of interacting particles [20,21,22]. In this manuscript we argue for a power law decrease similar to the Heisenberg limit for measurements of temperature with our interacting system.
The manuscript is organized as follows: In section 2 we outline the model and discuss the two limiting Ising and Dicke regimes and their phenomenology. In section 3 we review the phase diagram calculation [15] and we recast this calculation in terms of mean-field theory. In section 4 we discuss the phase diagram in detail. In Section V we analyze the application of the phase transition for metrology. Finally, we conclude in section 6.

The Ising Model
Consider a one-dimensional chain of N spins or two-level atoms realizing an Ising model with a transverse magnetic field, where h ≥ 0 is the transverse field and J ≥ 0 is the interaction strength between neighbouring spins and σ y N +1 ≡ σ y 1 . The model is not trivially easy to diagonalize, as the transverse field operators σ z i do not commute with the interaction terms σ y i σ y i+1 . In general, for h J the transverse field dominates in which case we expect σ y i ≈ 0 and σ z i = 1. For a strong coupling J h the interaction term favours parallel spin in the y-direction so σ z i ≈ 0 and the system possesses two states of equal energy with σ y i = +1, or σ y i = −1.
An (almost) exact diagonalization of the Ising Hamiltonian can be performed [23,24] using a Jordan-Wigner transformation which maps the spin-operators to fermionic operators. The particular mapping we employ here is c n = i j<n σ z j σ † n , where σ † n = (σ x n − iσ y n )/2 is the n'th spin-z step-up operator. Multiplying the spinoperators with phase factors, σ z j , depending on the spins at previous locations in the spin-chain leads to the fermionic anti-commutator relations c n , c † n = δ n,n .
The fermionic number operator is c † n c n = σ n σ † n = (1 − σ z n )/2 and the state with all spins pointing in the z-direction is mapped to the fermionic vacuum while c † n flips a spin from up to down thereby creating a fermion. The inverse mapping reads σ † n = −i j<n (1−2c † j c n )c n such that the y-y-interaction term in the Ising Hamiltonian becomes −Jσ y Due to the periodic boundary conditions we cannot remove the intermediate (1 − 2c † j c j )-factors in the final y-y-term −Jσ y N σ y 1 . This term can be written as and by adding and subtracting J(c † N − c N )(c † 1 + c 1 ) we can write the resulting Hamiltonian as where the last term is of relative order 1/N since 1 + j (1 − 2c † j c j ) is a projection operator onto the subspace of even total spin in the z-direction and therefore is of order 1.
By neglecting this less important final term the Hamiltonian is quadratic in the fermionic operators, and can be diagonalized using a Bogoliubov-transformation. The Bogoliubov-transformation can be decomposed into a single-particle basis change to momentum space and quasi-particle operators connecting particles of opposite momenta [23,24,25] and we get: where the operators γ k are the Bogoliubov transformed fermionic operators, (k) = 2(J 2 + h 2 − 2Jh cos(k)) 1/2 and k ∈ Z N has the form 2πn/N and is a reciprocal latticevector to the chain-lattice. The eigenmodes of H Ising correspond to free fermions with a dispersion relation given by (k). Note that (2) has a single unique ground state in contrast to (1) which has a two-fold degeneracy. This difference is due to the omitted term and will not be relevant in the thermodynamic limit.
Since the spectrum of the Ising model is so simple we can calculate the partition function for the system. The partition function is given by Z 0 Ising = Tr(exp(−βH Ising )) where β is the inverse temperature. We can easily calculate this trace, Using this result we can also calculate the free energy given by F Ising = −β −1 log Z 0 Ising . The diagonalization of H Ising is exact to order 1/N and to the same precision we can replace the sum over k ∈ Z N by an integral, Note that β (k, h, J)/2 = βJ(1+(h/J) 2 −2(h/J) cos(k)) 1/2 , and hence it is convenient to parametrize F Ising with the dimensionless quantitiesβ = βJ andh = h/J in the following. The Ising model exhibits an infinite order quantum phase-transition between a paramagnetic h J and ferromagnetic J h phase with critical point at J = h where a non-analyticity arises in (k) at k = 0. This system has been studied extensively [25] and is a prime example of a quantum phase transition. The quantum phase-transition does not survive to finite temperatures, but it has important consequences for the finite-temperature behaviour of the system near the critical point.

The Dicke model
Another phase transition model consists of spins or two-level atoms coupled to a harmonic oscillator mode with frequency ω, through the interaction where σ x i is the Pauli x-matrix acting on the i'th spin and a is the harmonic oscillator step down operator. We denote the coupling-strength of the oscillator to a single spin by g/ √ N . For a physical implementation with atoms inside a cavity with a mode volume V , the quantum field strength per photon is proportional to 1/ √ V , and our scaling thus corresponds to atoms with a constant spatial density which is well defined, also in the thermodynamic limit N → ∞.
The explicit form of the interaction is well established in quantum optical systems, and occurs both for two-level systems and for, e.g., Raman processes between two states via an intermediate excited state through absorption and stimulated emission of the quantum field and a classical control field. We note that in the Jaynes-Cummings model of a single two-level atom and a single field mode, the rotating wave approximation retains only the terms σ † a + σa † in V and provides a considerable simplification of the problem. The partition function for the Dicke-model has been calculated analytically in the thermodynamic limit and a second order phase transition has been identified [11]. While that calculation pertained to the rotating wave approximation it has been shown [12] that this does not change the Dicke phasetransition qualitatively. Within our application of a mean field approximation to the combined Dicke-Ising model we shall retain the full interaction (5) as, the rotating wave approximation is more difficult to deal with.
The Dicke model H osc + V can be realized experimentally as described in [14] where a dynamical version of the standard Dicke model is investigated in a cavity using four-level atoms coupled by Raman channels. The model parameters of this system are functions of the atomic and field parameters applied and can be tuned over large ranges.

The partition function for the Dicke-Ising model
The combined Dicke-Ising model has total Hamiltonian H = H Ising + H osc + V with three parameters g, h, and J which can be varied independently in atomic simulators of the model. With J = 0 the Hamiltonian realizes the Dicke model of N twolevel atoms interacting with a cavity-field via the dipole interaction. In this regime 2h correspond to the energy-splitting of the individual two-level atoms and g is the coupling strength to the cavity-field.

Coherent state integral
We will proceed by writing first the expression of the partition function for the full problem Z = Tr(exp(−βH)) with H = H Ising + H osc + V . Following [11] we write H as We will evaluate the trace over the oscillator-mode in the partition function using the coherent state representation of the field, and for this purpose we should bring exp(−βH) on normal ordered form with respect to the scaled operators The approximation of neglecting the commutator in the thermodynamic limit can be illuminated by Wick's theorem [26]: Any string of creation-and annihilation operators can be written as their normally ordered form plus additional terms. These terms are given by means of a contraction defined by C(AB) = AB−: AB:, where : AB: is the normal ordering of AB. Using this notation, Wick's theorem states that

Now, in our case
Hence, applying Wick's theorem to the expansion of the exponential in Z = Tr(exp(−βH)), and assuming that the limits in Z = lim N →∞ lim R→∞ R r=0 (−βH) r /r! can be interchanged, we see that all terms involving contractions will be of order 1/N or higher and can therefore be neglected in the thermodynamic limit. More generally any expression of the form Tr(f (b, b † )e −βH ), where f is a polynomial, can be calculated to an accuracy of 1/N using the normal order : f (b, b † )e −βH :. The validity of this truncation for the Dicke-model was discussed in [12].
When performing the trace of a normally ordered operator it is convenient to use the coherent states, |α , eigenstates of the annihilation operator: a|α = α|α , which yields The remaining trace over the spin degrees of freedom is exactly equivalent to the the original Ising model calculation, where the term involving the real part of the complex field argument acts as an additional magnetic field in the x direction, and hence the Ising model is biased by an effective magnetic field in the xz-plane with magnitude h 2 eff = h 2 + 4g 2 (α) 2 /N . Since this effective field lies within the plane orthogonal to the y-axis we can choose the direction of the effective field as a redefined z-axis and apply the same diagonalization as in the pure Ising problem.
This yields Z = d 2 α π exp(−βω |α| 2 + log Z 0 Ising (β, h eff , J)) and after performing the substitution z = α/ √ N we get where the dispersion relation (k, h eff (x), J) is for the effective magnetic field h eff and spin coupling J.

Mean field theory
We may also attack the original problem with an Ansatz replacing the interaction part V of the Hamiltonian by the mean field expression The last term compensates for double-counting of the interaction energy, while correlations in the mean-field fluctuations g With the Hamiltonian in this form the mean-fields split the Hamiltonian into two separate terms: a classically driven field mode H MF,osc = ωa † a + gN s x (a + a † )/ √ N and an Ising model with an additional magnetic field component along the x-direction, . This mean field Ising Hamiltonian has a field in the xz-plane with magnitude h 2 eff = h 2 + 4g 2 x 2 . The partition function for the mean field Hamiltonian therefore factors Z = Z MF,osc Z MF,Ising exp(2βgN s x x) and can be readily determined for arbitrary values of the mean field amplitudes x and s x . The corresponding free energy reads where f Ising = F Ising /N is the free energy per particle for the atoms as calculated in Eq.
(3). Note that the free energy includes a term representing the thermal distribution of the cavity-photons as well as the contribution depending on the field amplitude. The values of the mean-fields can now be obtained by minimizing the free energy. A short calculation reveals 1 N This shows that s x = −xω/g and that the Dicke-order parameter x should be found by minimizing ωx 2 + f Ising with respect to x. We will return to this minimization problem below.
In the mean-field description the interpretation of the physical properties of the system becomes clear and unambiguous. As an example, with the mean-field theory we can obtain expressions for various correlation functions for the atomic variables from the large amount of theory already present on the transverse Ising chain since in thermodynamic equilibrium the strongly coupled system is effectively identical to a rotated Ising chain.
It is reassuring, but hardly surprising, that the mean field result can be recovered from the coherent state integral for the thermodynamic limit N → ∞ result. By Laplace's (saddle point) method, one can replace the integral over coherent state amplitudes by discrete contributions from the location of the maximum of the exponential with respect to z in (6). It is easily shown that this maximization coincides with the mean field, identified by minimization of the free energy and by picking one of the maximizers in (6) to represent a symmetry-broken physical state of the system.

Phase diagram and analysis of the free energy
We now turn to the problem of finding the minimum of the free energy (7a),(7b) along with several important observables like the magnitude of the oscillator mean field, the spin magnetization and the susceptibility χ = ∂m z /∂h.
For fixed β and ω, the system is controlled by three parameters h, J and g. To illustrate the phase transitions in the system, we introduce a convenient way to plot different quantities as function of these variables in Fig. 1 and 2. In each plot, the sum h + J + g is fixed, and the corners of the triangles shown correspond to each of the three quantities acquiring the maximum value while the others vanish. The straight dotted lines converging to the corners of the triangles correspond to definite values of the ratio between the two quantities indicated on the edges of the triangles. These plots can be thought of as slices of the three-dimensional simplex defined by h + J + g + ω = /β where should be interpreted as the system energy-scale. In this coordinate representation, we show with colour coding the value of different interesting quantities.
Of particular interest in Fig. 1 is the oscillator field strength, represented by x 2 , and the susceptibility χ which are shown for the case where ω = 1, β = 100 and h + J + g = 1 for Fig. 1 and ω = 0.25, β = 4.0 and h + J + g = 1 for Fig. 2 where h, J and g are positive.
The edge of the triangle between g and h (i.e., with vanishing J) corresponds to the usual Dicke-model, while the edge between h and J (i.e., with vanishing g) corresponds to the usual Ising-model with the critical point at h/J = 1, showing up clearly as a signature in the variation of the susceptibility χ. We observe that this signature is present also for finite Dicke coupling parameter in the plot. The black curve in each plot shows where the Dicke phase transition occurs. In both Fig. 1 and Fig. 2 signatures of both first-and second-order phase transitions can be seen. Approximately below the line h/J = 1 the second order transition can be identified by the smooth increase in x 2 whereas above h/J = 1 one can discern a discontinuous jump in the order parameter.
The most significant difference between the cases presented in Figs. 1 and 2 is the susceptibility χ. For moderately low temperatures, β = 4 in Fig. 2, the signature of the Ising quantum phase transition is still clearly present, whereas for very low temperatures, β = 100 Fig. 1, the Ising phase transition becomes almost completely suppressed in the Dicke regime. Indeed, far into the Dicke regime (towards the right vertex in the triangles) the spin interactions do not appear to play any significant role.
Looking at the black and red solid lines in Fig. 1 one might be tempted to conclude that in the super-radiant phase we always have h eff /J > 1 (i.e. the red line does not penetrate into the area to the right of the black line). A close look at Fig. 2 will, however, reveal that this is not always the case. Indeed, for h + J + g = 0.6, ω = 0.62 and β = 4 there is a small part of the parameter-space where x > 0 and h eff < J. The physical reason why the Dicke phase transition almost, but not quite, suppresses the Ising transition remains to be understood. Formally, it occurs because the free energy has a minimum giving an effective magnetic field such that h eff /J > 1. The Ising critical point is therefore simply skipped in these cases and only an amputated signature of the Ising phase transition is present in cases such as shown in Fig. 1.
To investigate this phase transition in more detail, let us further consider the minimization of the free energy. The order-parameter x enters the mean field Isingterm via the effective magnetic field (h eff /J) 2 =h 2 + 4g 2 x 2 /J 2 . If we therefore introduce a rescaled order-parameterx = 2gx/J and a rescaled mode frequencỹ ω = ωJ/4g 2 we can write the free energy as a function of a few dimensionless quantities where C is a constant independent ofx. In order to understand the structure of the phase diagram it is necessary to investigate how the integral changes as a function ofβ = β/J andh compared toωx 2 . If we choose the variableh eff = h eff /J = h 2 + 4g 2 x 2 /J as the independent variable instead ofx, the integrand in (8) only depends onβ and the new variableh eff , whilẽ h eff ≥h imposes a boundary condition on the minimization with respect toh eff . To avoid confusion, we will consider F a function ofx and use the symbolF to denote the dependence onh eff . The system is in the super-radiant phase whenever the minimum inF occurs forh eff >h which implies x = 0. Examples ofF for representative values ofω andβ can be seen in figure 3. By a numerical investigation it is quickly revealed, thatF has at most a single local minimum (e.g. the curves (a), (b), (c) and (d) in figure 3) ath eff = 0 or no local minimum (curve (e) in figure 3). The existence and location of the minimum are thus solely determined byβ andω.
This implies that when keepingω andβ fixed the minimum of the free energy is either ath eff =h or at the local minimum ofF . If we imagine tuningh from high values towards low values (i.e, setting the boundary conditionh eff ≥h at different locations, for example along curve (c) in figure 3) the system will pass a second order phase-transition whenh passes the local minimum ofF . In the case that the local minimum is not the global minimum there will be ah > 0 whereF goes below its value at the local minimum implying that when one further lowersh the system will undergo a first-order phase transition into the normal state again. By the same reasoning ifF has a single global minimum (curve (a) and (b) in figure 3) there can only be a second order phase transition when tuningh.
The second-order phase transitions can be investigated in further detail using Gintzburg-Landau theory: In the neighbourhood of the second-order phase transition the order parameterx will always be small so we can expand the free energy F as a polynomial inx aroundx = 0: where I n is the n'th term in the Taylor expansion of the integral (8) with respect tox. The standard argument from Gintzburg-Landau theory is now that this fourth-order polynomial has a non-zero minimum whenω + I 2 (β,h) is negative. The second order phase-transition therefore occurs when I 2 (β,h) = −ω. By numerical investigation one finds that −I 2 is bounded by approximately 0.3356 implying that forω > 0.3356 no phase transition can occur. The first-order phase transition grows out of the second-order phase transition so there will be a region where the first-order jump in the order parameter is small. In that case we can still use Ginzburg-Landau theory and in particular we can find the point where the second-order transition changes to a first-order transition, i.e. when a local minimum in F changes from purely local to truly global. Again we analyze the polynomial expansion and one can show [27] that one needs to solve the system of equations I 4 (β,h) = 0 and I 2 (β,h) +ω = 0 to obtain the point where the phase transition changes nature. By investigating the functional form of I 4 it turns out that there is a minimal β c below which the first order phase-transition cannot occur. This value can be calculated numerically and is approximatelyβ c ≈ 1.1430. This identifies where the second-order transition changes to a first-order transition. To determine the first-order transition boundary for finite jumps in the order parameter, however, it is necessary to deal with the free energy F to all orders. Numerically it is not difficult to investigate for which value ofh the value ofF coincides with the value at the local minimum as described above. All this information has been combined into figure 4 where the phase-boundaries for various values ofω have been indicated. The coloured dotted lines represent first order transitions whereas the solid lines indicate second-order transitions. The black dashed curve indicates where I 4 (β,h) = 0 and its intercept with the curves I 2 (β,h) +ω = 0 indicates where the phase transition changes type between first-and second-order transitions.

Using a phase transition for a high precision measurement
A first order phase transition is interesting for many different reasons and here we consider its use as a measurement tool. Indeed, the standard description of a first order phase transition includes a discontinuous jump in the order parameter, and it is a relevant question, how precisely an experiment can locate the position of this discontinuous jump. The size of the jump-discontinuity usually scales linearly with the number of particles, while the width of the transition region often scales with an inverse power of this number, and under that assumption we shall present a simple model for the metrological sensitivity of the system. Since the phase transition occurs for rather non-trivial combinations of the temperature and the interaction parameters of the models, by changing some of these parameters in a controllable way, one may be able to select parameter ranges with particularly high sensitivity of the phase transition point to the value of the quantity being probed.
We will consider a measurement strategy where the system is probed at some range of values of some control parameter, e.g. a bias magnetic field. Fig. 5 shows how the cavity field order parameter varies as a function ofh for three different temperatures. To produce this plot, we have selected values ofω and thus an area of the phase diagram where the critical magnetic field depends strongly on the temperature, cf. the steepness of the dashed curves in Fig. 4. We expect that it will be possible to determine the critical value of the magnetic field with high precision, and since in this case a variation of the temperature of 1% changes the critical magnetic field by approximately 35%, the measurement of the critical field yields a very sensitive temperature measurement within the appropriate range of valuesβ ≈ 0.77 ± 1%. Sensitivity in, e.g., a lower temperature range is obtained if we chose a higher value ofω and scan a different range of values of the magnetic field.
The order parameter presented here is the intra cavity field intensity. We imagine  Fig. 4 where the slope of the phase transition line is large. For the blue curve 1/β = 0.77, for the green curve the temperature is 1% higher and for the red curve the temperature is 1% lower.
that the cavity leaks photons at a sufficiently low rate not to significantly disturb the thermodynamic steady state of the system, and herewith, detection of the intensity of the emitted light is a direct probe of the cavity field order parameter. We assume that the inverse temperatureβ is known to be close to some reference value, and we can then estimate the difference δβ by the best unbiased linear estimator as described in Appendix Appendix A. This estimator is given bŷ where n i is the detected number of photons in a given time while the controllable effective bias fieldh attains the valueh i . µ(h i ) is the expectation value of the photon number and σ 2 (h i ) is the photon number variance. µ (h) denotes the derivative of the expected photon number with respect to changes in inverse temperatureβ, and the expression applies within a narrow range where a linear variation of the expected photon number withβ is valid.
In the limit of high bias field resolution, the sum in the estimator can be converted into an integral, and one can determine the variance of the estimate: Var(δβ(n(h)) = 1/ µ (h) 2 /σ 2 (h)dh, see details in the appendix.
So far the arguments have been of a general nature. Let us now assume the Dicke-Ising model, in which the photon number distribution is well described as a thermal state below and a displaced thermal state above the Dicke phase transition. The first and second moments of such distributions can be calculated using, e.g. the positive P-representation for the thermal state, µ(h,β) = n = N x 2 +n σ 2 (h,β) = Varn = N x 2 (1 + 2n) +n +n 2 Recall that the order parameter x 2 is a function of the system parametersh,β and ω and in the thermodynamic limit it has a discontinuous jump at the dashed lines shown in Fig. 4. For a finite system, however, the phase transition constitutes a smooth curve with a fast increase of the order parameter. The width of this region is not easy to determine but finite size effects in phase transitions tend to smoothen phase transitions leading to a decreasing width as N increases. Indeed, [28] finds that, in general, the width scales as 1/N , but for the sake of generality we assume a scaling N −γ , γ > 0.
Since both n and Varn are proportional to N x 2 both µ(h), µ (h) and σ 2 (h) will carry a signature of this power law. To be explicit assume that Nx 2 (h,β) = N tanh(N γ (h −h c (β)) whereh c is the critical value ofh c as a function ofβ. Then The function µ only has support in a region of width 1/N γ nearh 0 ≡h c (β 0 ). The variance of the estimatê This is our main result of this section, showing that the sensitivity is better than the "standard limit" where one expects Varβ ∼ 1/N , and depending on the character of the finite size effects (the power γ), it is potentially also better than the Heisenberg detection limit. With the result from [28], γ = 1 the accuracy is actually at the Heisenberg limit. Note that the above argument is quite general and applies to any first order transition with an intensive order parameter.
The termh c is included in our expression in order to show explicitly that the sensitivity depends on the curve of critical points in the phase diagram of the system. From Fig. 4 we see thath c can be chosen large for arbitrarily small temperatures by tuningω. The large value ofh c , however, comes at a cost: The slope of h c is highest near the thick dashed curve, which is also where the first order transition has small amplitude and changes to a second order transition. In a concrete implementation, the values ofω and the range of effective magnetic fields need to be chosen with care to reflect the actual scaling N −γ and the size of the jump discontinuity.
With an adaptive measurement scheme, we imagine that the number of iterations with differenth i for a reliable detection of the critical value of theh-parameter can be optimized. It is clear that a more detailed investigation is necessary in order to quantify the accuracy and scaling of resources of such measurements. Indeed, the specific power law 1/N γ for the transition width is only a convenient Ansatz, and a non-mean field calculation on a finite system will be needed in order to investigate the approach towards the thermodynamic limit in more detail. Furthermore, the critical properties and the long range correlations of the system may possibly lead to even better estimates by use of the recent techniques of quantum non-linear parameter estimation [29,21]. These issues we shall defer to a later publication.

Conclusion
We have investigated the thermodynamic properties of a Dicke-Ising model incorporating both the quantum transverse Ising model and the Dicke model as special limiting cases. We have derived expressions for the free energy using a coherent state integral similar to [11] but also using a mean field theory with a clearer interpretation for the field statistics. The combined model exhibits a first order phase transition which is not present in either of the two separate models. By a simple numerical search the free energy minimum can be identified and the value of the Dicke mean field and the magnetic susceptibility can be determined as functions of all physical parameters of the model, cf. Figs. 1 and 2.
Using the free energy and Ginzburg-Landau theory we also investigated the complete phase diagram as shown in Fig. 4. The Dicke phase transition occurs also for moderate inter-particle interactions and the Ising phase transition is also well preserved for weak and moderate light-matter couplings. In a small area of the parameter space both phase transitions coexist closely together, but for a stronger Dicke model interaction the resulting mean field puts the system in a regime without any observable Ising phase transition.
The Dicke-Ising model constitutes an interesting mix of second-, first-and infinite order phase transitions. The interplay of these phase transitions and a complete description beyond the mean field approximation of the fundamental excitations at the critical points would be and interesting continuation along the lines of this work. In addition to its fundamental theoretical interest, the first order phase transitions provides a tool for precise measurements of, e.g. the magnetic bias field or of the temperature. We have presented a simple estimate of the accuracy of such a measurement device showing that the variance scales as 1/N 1+γ which is better than the standard limit 1/N for independent measurements on N particles. This work was supported by the European Union Integrated Project AQUTE.

Appendix A. Best unbiased linear estimator
We will consider the best unbiased linear estimator in a situation where an experimenter performs a sequence of measurements where she scans a parameter x (e.g. a bias magnetic field) in order to uncover another, unknown, parameter q (e.g. the location of a critical point). The experimenter measures a discrete stochastic variable n (e.g. photon number) which has probability distribution p(n; x, q). In the measurement the experimenter thus collects, for a fixed q, the values n i corresponding to selected x i . In the following we assume that the unknown parameter q is close to a reference value, which, without loss of generality, we take to be zero. For q 1 we can then expand the moments of p in a Taylor expansion such that where we have expanded the mean µ(x, q) to first order (and µ(x) ≡ µ(x, q), µ (x) ≡ (∂ q µ)(x, 0)) and the variance σ 2 (x, q) to zeroth order in q.
A linear estimator is of the form where {g i } are weighting coefficients, {n i } are the observed values of n i corresponding to the chosen values x i and c is a constant. To find an unbiased estimator we require E [q({n i })] ∝ q which implies c = − i g i µ 0 (x i ). To find the best linear estimator we optimize the signal-to-noise ratio E [q({n i })] 2 /Var(q({n i })) with respect to the vector g i . To second order in q the signal-to-noise ratio is q 2 ( i g i µ (x i )) 2 / i g 2 i σ 2 (x i ). The minimum of the signal-to-noise ratio is obtained when g i ∝ µ (x i )/σ 2 (x i ).
The constant of proportionality should then be chosen such that q = E [q({n i })] = i g i (E [n i ] − µ(x i )) which gives the condition i g i µ (x i ) = 1 and the normalization constant is given by A −1 = i (µ i ) 2 /σ 2 i . With this normalization the variance of the estimate q is