/******************************************************************************* * This file is part of SWIFT. * Copyright (c) 2019 Joel Pfeffer (j.l.pfeffer@ljmu.ac.uk) * * This program is free software: you can redistribute it and/or modify * it under the terms of the GNU Lesser General Public License as published * by the Free Software Foundation, either version 3 of the License, or * (at your option) any later version. * * This program is distributed in the hope that it will be useful, * but WITHOUT ANY WARRANTY; without even the implied warranty of * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the * GNU General Public License for more details. * * You should have received a copy of the GNU Lesser General Public License * along with this program. If not, see . * ******************************************************************************/ /* Standard includes */ #include /* Standard includes */ #include #include /* Local includes */ #include "dsyevj3.h" #include "engine.h" #include "stars_io.h" /** * @brief Free-fall time as a function of density * * @param rho The gas density in kg/m^3 * @param G Gravitational constant * @return Free-fall time */ __attribute__((const)) INLINE static double f_tff(const double rho, const double G) { /* free-fall time sqrt(3 * pi / (32 * G * rho) */ return sqrt(0.09375 * M_PI / (G * rho)); } /** * @brief Compute the Mach number * * @param sigmaloc Velocity dispersion * @param csloc Sound speed * @return Mach number */ __attribute__((const)) INLINE static double f_mach(const double sigmaloc, const double csloc) { return sigmaloc / csloc; } /** * @brief Dispersion of overdensity PDF as a function of Mach number and * magnetic pressure ratio, based on e.g. Padoan & Nordlund (2011). * * @param mach Mach number * @param beta0 Magnetic pressure ratio * @return Dispersion of overdensity PDF */ __attribute__((const)) INLINE static double f_sigrho(const double mach, const double beta0) { /* Constant, e.g. Padoan & Nordlund 2002 */ const double b = 0.5; return sqrt(log(1. + 3. * b * b * mach * mach * beta0 / (beta0 + 1.))); } /** * @brief Critical overdensity for star formation in the sSFR_ff. * * Equation 27 of Krumholz & McKee, 2015, ApJ, 630, 1, pp. 250-268. * * @param qvir GMC virial ratio * @param mach Mach number * @return Crititcal overdensity */ __attribute__((const)) INLINE static double f_xcrit(const double qvir, const double mach) { /* constant (table 2) */ const double phix = 1.12; return (1. / 15.) * M_PI * M_PI * phix * phix * qvir * mach * mach; } /** * @brief Specific star formation rate per free-fall time. * * Uses Eq. 20 of Krumholz & McKee, 2015, ApJ, 630, 1, pp. 250-268 * when sflaw == 1 * OR * Eq. 13 of Elmegreen, 2002, ApJ, 577, 1, pp. 206-220 * which is just a constant value. * * @param mach Mach number * @param props Properties of the stars model * @return Specific star formation rate per free-fall time */ __attribute__((const)) INLINE static double f_sfrff( const double mach, const struct stars_props* props) { /* Star formation law */ if (props->sflaw == 1) { /* Krumholz & McKee (2005). Depends on parameters: * qvir - GMC virial ratio * beta0 - Magnetic pressure ratio * ecore - Formation efficiency of protostellar cores */ /* constant (table 2) */ const double phit = 1.91; /* critical overdensity for star formation */ const double xcrit = f_xcrit(props->qvir, mach); /* overdensity PDF dispersion */ const double sigrho = f_sigrho(mach, props->beta0); const double erf_term = (-2. * log(xcrit) + sigrho * sigrho) / (2. * M_SQRT2 * sigrho); return 0.5 * props->ecore * (1. + erf(erf_term)) / phit; } else { /* Elmegreen (2002) */ return props->const_sfrff; } } /** * @brief Overdensity PDF of the ISM * * @param x Overdensity * @param mulnx Logarithmic mean * @param sig Dispersion * @return Overdensity PDF */ __attribute__((const)) INLINE static double f_dpdx(const double x, const double mulnx, const double sig) { const double dif = log(x) - mulnx; /* 1. / (sqrt(2 pi) * sig * x) * exp (-dif^2 / 2 sig^2) * * Note: 1 / sqrt(2 pi) = (1/4) * (2 / sqrt(pi)) * sqrt(2) */ return ((0.25 * M_2_SQRTPI * M_SQRT2) / (sig * x)) * exp(-0.5 * dif * dif / (sig * sig)); } /** * @brief Estimate for gas surface density, assuming an equilibrium disc * * @param rholoc Gas volume density * @param sigmaloc Gas velocity dispersion * @param fgas Local gas fraction * @param G Gravitational constant * @return Gas surface density */ __attribute__((const)) INLINE static double f_surfg(const double rholoc, const double sigmaloc, const double fgas, const double G) { return sqrt(2. * fgas * rholoc * sigmaloc * sigmaloc / (M_PI * G)); } /** * @brief Calculate local bound fraction of star formation * * @param epsilon Gas volume density * @param logZ log10 of stellar metal mass fraction Z * @param props Properties of the stars model * @return local bound fraction */ __attribute__((const)) INLINE static double f_fbound( const double epsilon, const double logZ, const struct stars_props* props) { double bound; if (props->CFE_met_dep_fbound) { /* Metallicity dependence for fbound-SFE, based on Grudic et al. 2021. * Approximate fit to simulations in their figure 3 with a double power * law (with indices alpha and 0) and shape parameter epsilon^beta. * "Saturation" value from Li et al. 2019 to shift vertical scale. * Linearly interpolate the power-law index in log(Z) so that * fbound ~ Z^beta (with beta a function of epsilon), consistent with a * power-law scaling for wind mass loss rates (e.g. Vink et al. 2001). The * lower metallicity value sets the lower alpha limit based on the * "no winds" run in Grudic et al. 2021. */ double alpha; if (logZ < props->fbound_logZ[0]) { /* Lower metallicity value sets lower limit on alpha */ alpha = props->fbound_alpha[0]; } else { /* Interpolate the power law index by particle metallicity */ alpha = props->fbound_alpha[0] + (logZ - props->fbound_logZ[0]) * (props->fbound_alpha[1] - props->fbound_alpha[0]) / (props->fbound_logZ[1] - props->fbound_logZ[0]); } bound = props->fbound_saturation * pow((1. + props->fbound_scale) / (1. + props->fbound_scale / pow(epsilon, props->fbound_beta)), alpha); } else { /* Original one-to-one assumption in K12 */ bound = epsilon; } return bound; } /** * @brief Naturally bound fraction of star formation * * @param rholoc Gas volume density * @param sigmaloc Gas velocity dispersion * @param mach Mach number * @param xloc Overdensity * @param surfGMC Giant molecular cloud surface density * @param sfrff Specific star formation rate per free-fall time * @param props Properties of the stars model * @param phys_const Physical constants * @return Local SFE */ static double f_fstar(const double rholoc, const double sigmaloc, const double mach, const double xloc, const double surfGMC, const double sfrff, const struct stars_props* props, const struct phys_const* phys_const) { /* Extract some physics constants (in internal units) */ const double k_B = phys_const->const_stefan_boltzmann; const double c = phys_const->const_speed_light_c; const double G = phys_const->const_newton_G; /* Extract some model constants */ const double t_SN = props->tsn; const double t_view = props->tview; const double phi_fb = props->const_phifb; const double kappa0 = props->const_kappa0; const double psi = props->const_psi; const double ecore = props->ecore; /* trapping ratio */ const double phitrap = 0.2; /* Free-fall time */ const double t_ff = f_tff(xloc * rholoc, phys_const->const_newton_G); /* Star Formation Efficiencies */ /* SN feedback */ double efb; if (props->radfb == 0 || props->radfb == 2) { efb = 0.5 * sfrff * (t_SN / t_ff) * (1. + sqrt(1. + 4. * t_ff * sigmaloc * sigmaloc / (phi_fb * sfrff * t_SN * t_SN * xloc))); } else efb = 1.; /* radiative feedback */ double efbrad; if (props->radfb > 0.) { /* surface density on which radiative feedback acts */ const double surffb3 = surfGMC * surfGMC * surfGMC; const double factor = phitrap * kappa0 * kappa0 / k_B; efbrad = 2. / (factor * psi * surffb3) * (sqrt(1. + 2. * M_PI * c * G * factor * surfGMC * surffb3) - 1.); } else efbrad = 1.; /* star formation is incomplete/still ongoing */ const double einc = sfrff * t_view / t_ff; /* Local SFE is minimum of all efficiencies: * min[maximum, SNfeedback, radiative feedback, incomplete] */ return min4(ecore, efb, efbrad, einc); } /** * @brief Calculate the fraction of star formation occurring in bound clusters * i.e. the "cluster formation efficiency" * Kruijssen, J. M. D., 2012, MNRAS 426, 3008 * * @param rholoc Local gas volume density in kg/m^3 * @param sigmaloc Local 1D gas velocity dispersion in m/s * @param csloc Local gas sound speed in m/s * @param fgas Local gas fraction * @param logZ log10 of stellar metal mass fraction Z * @param props Properties of the stars model * @param phys_const Physical constants * @return The CFE */ static double f_cfelocal(const double rholoc, const double sigmaloc, const double csloc, const double fgas, const double logZ, const struct stars_props* props, const struct phys_const* phys_const) { /* Parameters: * sflaw - star formation law: 0=Elmegreen (2002, default) * 1=Krumholz & McKee (2005) * qvir - giant molecular cloud virial parameter (default=1.3) * tsn - time of the first supernova (default=3.*Myr) * tview - time at which CFE is determined (default=10.*Myr) * surfGMC - giant molecular cloud surface density (default=100.*Msun/pc^2.) * ecore - maximum (protostellar core) star formation efficiency (default=.5) * beta0 - turbulent-to-magnetic pressure ratio (default=1.e10) * radfb - feedback: 0=supernova, 1=radiative, 2=both (default=0) * cce - Include cruel cradle effect? (default=0) * xsurv - Critical overdensity for cruel cradle effect (default=0.) */ /* Extract some physics constants (in internal units) */ const double G = phys_const->const_newton_G; /* Integrate overdensity PDFs to obtain bound fraction of star formation. * Number of integration steps (checked to be sufficient for convergence) */ const int nx = 1000; /* estimate of gas surface density */ const double surfg = f_surfg(rholoc, sigmaloc, fgas, G); const double surfGMC = (surfg > props->surfGMC_min) ? surfg : props->surfGMC_min; /* mach number */ const double mach = f_mach(sigmaloc, csloc); /* specific star formation rate per free-fall time */ const double sfrff = f_sfrff(mach, props); /* dispersion of overdensity PDF */ const double sigrho = f_sigrho(mach, props->beta0); /* logarithmic mean of overdensity PDF */ const double mulnx = -0.5 * sigrho * sigrho; /* Minimum overdensity. If calculating the cruel cradle effect and critical * overdensity below minimum, then adjust */ double xmin1 = exp(mulnx - 5. * sigrho); if (props->cruel_cradle > 0 && props->xsurv < xmin1) xmin1 = props->xsurv; /* Maximum overdensity. If calculating the cruel cradle effect and critical * overdensity below maximum, then adjust */ double xmax = exp(mulnx + 10. * sigrho); if (props->cruel_cradle > 0 && props->xsurv > xmax) xmax = props->xsurv; /* initialize first integral */ const double ratio1 = xmax / xmin1; const double step_size_term1 = pow(ratio1, 1. / (2. * nx)) - pow(ratio1, -1. / (2. * nx)); double f1 = 0.; /* denominator integral */ for (int ix = 0; ix < nx; ix++) { /* Overdensity */ const double xloc = xmin1 * pow(ratio1, (ix - 0.5) / nx); /* Step size */ const double dx = xloc * step_size_term1; /* Local SFE */ const double fstar = f_fstar(rholoc, sigmaloc, mach, xloc, surfGMC, sfrff, props, phys_const); /* Local bound fraction */ double bound; if (props->cruel_cradle == 0 || props->cruel_cradle == 2) { /* If not calculating the cruel cradle effect but the naturally bound * fraction of SF, the denominator should contain all SF. * If calculating the cruel cradle effect with respect to all SF, the * denominator should contain all SF. */ bound = 1.; } else bound = f_fbound(fstar / props->ecore, logZ, props); /* integral part 1 */ const double integral = bound * fstar * xloc; /* overdensity PDF, i.e. integral part 2 */ const double dpdx = f_dpdx(xloc, mulnx, sigrho); /* integral */ f1 += integral * dpdx * dx; } /* if calculating the cruel cradle effect set minimum overdensity to critical * overdensity */ const double xmin2 = (props->cruel_cradle > 0) ? props->xsurv : xmin1; /* initialize second integral */ const double ratio2 = xmax / xmin2; const double step_size_term2 = pow(ratio2, 1. / (2. * nx)) - pow(ratio2, -1. / (2. * nx)); double f2 = 0.; /* numerator integral */ for (int ix = 0; ix < nx; ix++) { /* overdensity */ const double xloc = xmin2 * pow(ratio2, (ix - 0.5) / nx); /* Step size */ const double dx = xloc * step_size_term2; /* local SFE */ const double fstar = f_fstar(rholoc, sigmaloc, mach, xloc, surfGMC, sfrff, props, phys_const); /* local bound fraction */ double bound; if (props->cruel_cradle == 2) { /* if calculating the cruel cradle effect with respect to all SF, the * numerator should contain all SF */ bound = 1.; } else bound = f_fbound(fstar / props->ecore, logZ, props); /* integral part 1 */ const double integral = bound * fstar * xloc; /* overdensity PDF, i.e. integral part 2 */ const double dpdx = f_dpdx(xloc, mulnx, sigrho); /* integral */ f2 += integral * dpdx * dx; } /* The cluster formation efficiency */ if (f1 == 0) return 0.; else return f2 / f1; } /** * @brief Local supernovae feedback timescale * * @param rholoc Local gas volume density * @param sigmaloc Local 1D gas velocity dispersion * @param csloc Local gas sound speed * @param props Properties of the stars model * @param phys_const Physical constants * @return The supernovae feedback timescale of molecular cloud */ static double feedback_timescale(const double rholoc, const double sigmaloc, const double csloc, const struct stars_props* props, const struct phys_const* phys_const) { /* Overdensity */ /* const double xloc = 1.; */ /* Mach number */ const double mach = f_mach(sigmaloc, csloc); /* specific star formation rate per free-fall time */ const double sfrff = f_sfrff(mach, props); /* free-fall time */ const double tff = f_tff(rholoc, phys_const->const_newton_G); const double tsn = props->tsn; const double phi_fb = props->const_phifb; /* SN feedback timescale */ return 0.5 * tsn * (1. + sqrt(1. + 4. * tff * sigmaloc * sigmaloc / (phi_fb * sfrff * tsn * tsn))); } /** * @brief Swap two numbers. */ __attribute__((always_inline)) INLINE static void swap(double* a, double* b) { double tmp = *a; *a = *b; *b = tmp; } /** * @brief Sort array with length 3 * * @param val Length 3 array to sort */ __attribute__((always_inline)) INLINE static void sort3(double val[3]) { if (val[0] > val[1]) swap(&val[0], &val[1]); if (val[0] > val[2]) swap(&val[0], &val[2]); if (val[1] > val[2]) swap(&val[1], &val[2]); } /** * @brief Local cluster evaporation timescale * * @param sp The #spart to act upon. * @param props Properties of the stars model. * @param a3_inv Inverse cube of the current expansion factor * @return Cluster evaporation timescale */ double get_local_evaporation_t0(struct spart* restrict sp, const struct stars_props* props, const double a3_inv) { /* Convert tensor to physical units. Note the lower elements are never * referenced by dsyevj3 */ double tensor[3][3]; tensor[0][0] = (double)sp->tidal_tensor[2][0] * a3_inv; tensor[1][1] = (double)sp->tidal_tensor[2][1] * a3_inv; tensor[2][2] = (double)sp->tidal_tensor[2][2] * a3_inv; tensor[0][1] = (double)sp->tidal_tensor[2][3] * a3_inv; tensor[0][2] = (double)sp->tidal_tensor[2][4] * a3_inv; tensor[1][2] = (double)sp->tidal_tensor[2][5] * a3_inv; /* get the eigenvalues, then sort */ double tide_val[3]; dsyevj3(tensor, tide_val); sort3(tide_val); double Omega2; if (props->Omega_is_lambda2 && tide_val[2] < 0) { /* Correct version of Omega for a spherical distribution, in principle. * When tideval[2] > 0 we can't guarantee that tideval[1] will be < 0 */ Omega2 = fabs(-tide_val[1]); } else { /* Version in E-MOSAICS */ Omega2 = fabs(-tide_val[0] - tide_val[1] - tide_val[2]) * (1. / 3.); } /* Strength of the tidal field, T + Omega^2 */ const double tidal_strength = tide_val[2] + Omega2; /* Compute imescale of evaporation in the tidal field */ if (tidal_strength > 0.) { /* Interpolate between King W0=5 and W0=7 profiles */ const double t0evapsun = props->const_t0evapsunW5 + 0.5 * (props->W0 - 5.) * (props->const_t0evapsunW7 - props->const_t0evapsunW5); return t0evapsun / sqrt(tidal_strength / props->const_tidesun); } else { /* Compressive field */ return DBL_MAX; } } /** * @brief Calculate the tidal heating parameters * * @param sp The particle to act upon. * @param age Age of star. * @param dt Particle timestep. * @param time Current time. * @param a3_inv Inverse cube of the current expansion factor * @param Mlheatsum (return) Tidal heating for shocks at this timestep. */ void calc_heating(struct spart* restrict sp, const double age, const double dt, const double time, const double a3_inv, double* Mlheatsum) { /* fraction of maximum at which duration is measured */ const double sigma1_width = 0.88; const double sigma2_width = 0.61; /* Time at the previous step */ const float tlast = (float)(time - dt); /* loop over tidal tensor components, to determine minima, maxima, tshock * and tidal heating for each component. * Note we only store upper six components of tidal tensor. */ for (int i = 0; i < 6; i++) { const double tide1 = (double)sp->tidal_tensor[0][i] * a3_inv; const double tide2 = (double)sp->tidal_tensor[1][i] * a3_inv; const double tide3 = (double)sp->tidal_tensor[2][i] * a3_inv; int apply_shock = 0; /* Is the first point for this particle is a maximum? */ if (age == dt && fabs(tide3) < fabs(tide2)) { /* Time of last maximum is previous time step */ sp->tmaxsh[i] = tlast; sp->tidmax[i] = tide2; sp->extreme_max[i] = 1; sp->shock_duration[i] = 0.; sp->shock_indicator[i] = 0; } /* If previous time step marks a valid minimum */ if (fabs(tide2) < fabs(tide1) && fabs(tide2) < fabs(tide3) && fabs(tide2) < 0.88 * fabs(sp->tidmax[i])) { /* if previous extreme was a maximum, or this minimum is lower than * previous */ if (sp->extreme_max[i] || (fabs(tide2) < fabs(sp->tidmin[i]))) { /* no shock duration is present (shallow max) */ if (sp->shock_indicator[i] == 0) sp->shock_duration[i] = 0.5 * (tlast - sp->tminsh[i]); /* Time of last minimum was previous time step */ apply_shock = 1; sp->tminsh[i] = tlast; sp->tidmin[i] = tide2; sp->extreme_max[i] = 0; } } /* If previous time step marks a valid maximum */ if (fabs(tide2) > fabs(tide1) && fabs(tide2) > fabs(tide3) && 0.88 * fabs(tide2) > fabs(sp->tidmin[i])) { // if previous extreme was a minimum, or this maximum is higher than // previous if (!sp->extreme_max[i] || (fabs(tide2) > fabs(sp->tidmax[i]))) { /* Time of last maximum was previous time step */ sp->tmaxsh[i] = tlast; sp->tidmax[i] = tide2; sp->extreme_max[i] = 1; sp->shock_duration[i] = 0.f; sp->shock_indicator[i] = 0; } } /* First order shock duration */ if (sp->extreme_max[i] && sp->shock_indicator[i] == 0 && fabs(tide3) < sigma1_width * fabs(sp->tidmax[i])) { sp->shock_indicator[i] = 1; /* time interpolation fraction */ const double dt_fraction = (sigma1_width * fabs(sp->tidmax[i]) - fabs(tide2)) / (fabs(tide3) - fabs(tide2)); /* interpolate to time of shock duration measurement */ const double shock_time = time + (dt_fraction - 1.) * dt; /* set shock duration indicator to total width of 1 sigma */ /* 2 * dt / sigma(=1) */ sp->shock_duration[i] = 2. * (shock_time - sp->tmaxsh[i]); } /* Second order shock duration */ if (sp->extreme_max[i] && sp->shock_indicator[i] < 2 && fabs(tide3) < sigma2_width * fabs(sp->tidmax[i])) { sp->shock_indicator[i] = 2; /* time interpolation fraction */ const double dt_fraction = (sigma2_width * fabs(sp->tidmax[i]) - fabs(tide2)) / (fabs(tide3) - fabs(tide2)); /* interpolate to time of shock duration measurement */ const double shock_time = time + (dt_fraction - 1.) * dt; /* set shock duration indicator to total width of 2 sigma */ /* 2 * dt / sigma(=2) */ sp->shock_duration[i] = shock_time - sp->tmaxsh[i]; } /* If a shock was completed during the previous time step and mass loss * should be applied */ if (apply_shock) { /* Tidal heating for the soon-to-occur application of mass loss * Note this could be negative, but is later applied as heatsum^2 */ Mlheatsum[i] = sp->heatsum[i]; /* reset and start integrating heating for next shock */ sp->shock_indicator[i] = 0; sp->heatsum[i] = 0.5 * (tide2 + tide3) * dt; } else { /* Keep summing tidal heating */ sp->heatsum[i] += 0.5 * (tide2 + tide3) * dt; } } /* For each tensor component */ } /** * @brief Compute the Weinberg adiabatic correction * * @param mass Cluster mass in Solar masses. * @param rh Cluster half-mass radius in pc. * @param tshock Tidal shock duration. * @param W0 King profile parameter. * @param etapot Potential energy constant. * @param G Gravitational constant. * @return Adiabatic correction. */ double adiabatic_correction(const double mass, const double rh, const double tshock, const double W0, const double etapot, const double G) { /* constant in omega for King profile W0=[5,7] */ const double cst5 = 0.68; const double cst7 = 0.82; /* Interpolate */ const double cst = cst5 + 0.5 * (W0 - 5.) * (cst7 - cst5); /* angular frequency of stars */ const double omega = cst * sqrt(8. * etapot * etapot * etapot * G * mass / (rh * rh * rh)); /* number of stellar revolutions per passage */ const double x = omega * tshock; const double Aw = 1. / ((1. + x * x) * sqrt(1. + x * x)); #ifdef SWIFT_DEBUG_CHECKS if (!isfinite(Aw)) error("Adiabatic correction is not finite!"); #endif return Aw; } /** * @brief Compute star cluster mass loss by tidal shocks. * * @param sp The #spart. * @param props Properties of the stars model. * @param Mlheatsum Tidal heating to apply at this timestep. * @param mass Star cluster mass. * @param G Gravitational constant. * @return Mass loss by tidal shocks. */ double mass_loss_shocks(struct spart* restrict sp, const struct stars_props* props, const double* Mlheatsum, const double mass, const double G) { /* Anything to do? */ double max_heatsum = 0.; for (int i = 0; i < 6; i++) { if (fabs(Mlheatsum[i]) > max_heatsum) max_heatsum = fabs(Mlheatsum[i]); } if (max_heatsum == 0.) return 0.; double max_shock_duration = 0.; for (int i = 0; i < 6; i++) { if (fabs(sp->shock_duration[i]) > max_shock_duration) max_shock_duration = fabs(sp->shock_duration[i]); } if (max_shock_duration == 0.) return 0.; /* Ok, now complete the shock */ /* potential energy constant */ const double etapot = 0.4; /* fraction of shock energy input used in mass loss */ const double efrac = 0.25; /* mean(M(r)*r)/(Mtot*rh) for King profile W0=[5,7] */ const double zeta5 = 0.81034; const double zeta7 = 1.0273; /* ratio of square radii for King profile W0=[5,7] */ const double r2avrh25 = 2.; const double r2avrh27 = 3.5; /* Now interpolate */ const double zeta = zeta5 + 0.5 * (props->W0 - 5.) * (zeta7 - zeta5); const double r2avrh2 = r2avrh25 + 0.5 * (props->W0 - 5.) * (r2avrh27 - r2avrh25); /* Half-mass radius */ const double rh = props->rh; /* Sum up the tidal heating */ double heatparam = 0.; for (int i = 0; i < 6; i++) { /* Weinberg adiabatic correction */ double Aw = adiabatic_correction(mass, rh, (double)sp->shock_duration[i], props->W0, etapot, G); /* Add the mirror for off-diagonal tensor terms */ double factor = 1.; if (i > 2) factor = 2.; heatparam += factor * Mlheatsum[i] * Mlheatsum[i] * Aw; } /* first order tsh/P */ const double t1 = 3. * etapot * G * mass / (efrac * rh * rh * rh * r2avrh2 * heatparam); /* factor to include second order energy gain */ const double InclE2 = 1. + (12. * zeta) / (5. * etapot * r2avrh2); /* Disruption time / shock period: tsh_int = t1 / InclE2 */ return -mass * InclE2 / t1; } /** * @brief Compute star cluster mass loss by evaporation. * * @param props Properties of the stars model. * @param mass Star cluster mass. * @param initial_mass Initial star cluster mass. * @param t0 Evaporation timescale. * @param dt Particle timestep. * @param logZ log10 of stellar metal mass fraction Z. * @param G Gravitational constant. * @return Mass loss by evaporation. */ double mass_loss_evaporation(const struct stars_props* props, const double mass, const double initial_mass, const double t0, const double dt, const double logZ, const double G, const double const_solar_mass) { /* Dissolution law exponent for W0 = [5,7] King profile */ const double gamma5 = props->gamma_W5; const double gamma7 = props->gamma_W7; /* adapt dissolution law exponent to King profile parameter */ const double gamma_min = gamma5 + 0.5 * (props->W0 - 5.) * (gamma7 - gamma5); double gamma = gamma_min; if (props->met_dep_evaporation) { /* Add dependence on metallicity, following Gieles & Gnedin 2023 */ const double gamma_max = props->gamma_max_Zdep; const double logZ_ymax = props->logZ_gamma_max; const double logZ_ymin = props->logZ_gamma_min; gamma = gamma_min + (logZ - logZ_ymin) * (gamma_max - gamma_min) / (logZ_ymax - logZ_ymin); if (gamma < gamma_min) gamma = gamma_min; else if (gamma > gamma_max) gamma = gamma_max; } double dm = 0.; /* Add the 'isolated' evaporation term, following Gieles & Baumgardt 08 */ if (props->spitzer_evap_term) { /* Mass loss fraction per relaxation time in 'isolated' regime * (Spitzer 87, Gieles & Baumgardt 08) */ const double xi0 = 0.0074; /* Mean star mass for Kroupa IMF. Convert to internal units */ const double mmean = 0.64 * const_solar_mass; /* Cluster half-mass radius */ const double rh = props->rh; /* Relaxation time (Spitzer & Hart 71; Gierz & Heggie 94) */ const double N = mass / mmean; const double trh = 0.138 * rh * sqrt(rh * N / (mmean * G)) / log(0.11 * N); dm += -xi0 * mass * dt / trh; } /* Fraction of evaporation timescale. Includes internal unit correction */ double dt_t0; if (gamma != gamma_min) { /* Gieles & Gnedin 2023 rewrite * Mdot \propto -Mconst * (M/Mconst)^(1-y) * as * Mdot \propto -Mconst * (M/Minit)^(1-y) * (Minit/Mconst)^(1-x) * where x = 2/3 */ dt_t0 = dt / t0 * pow(const_solar_mass, gamma_min) * pow(initial_mass, gamma - gamma_min); } else { dt_t0 = dt / t0 * pow(const_solar_mass, gamma); } /* First and second order terms */ double dm_order1 = -pow(mass, 1. - gamma) * dt_t0; double dm_order2 = (1. - gamma) * pow(mass, 1. - 2. * gamma) * dt_t0 * dt_t0; /* Was it mass loss? Set to first order if mass gain. * dm>0 only occurs just before disruption because we don't consider third * order (etc) effects */ if (dm_order1 + dm_order2 > 0.) dm += dm_order1; else dm += dm_order1 + dm_order2; return dm; } /** * @brief Normalisation of Schechter function for probability function */ double schechternorm(double mmin, double mmax, double slope, double mstar) { /* integration resolution */ const int nintsteps = 10000; /* starting mass */ double mc = mmin; double norm = 0.; /* Constant integration step in log space */ const double dm = pow(mmax / mmin, 1. / (double)nintsteps) - 1.; /* integrate to determine norm (total number) */ for (int n = 0; n < nintsteps; n++) { /* dm */ double dmc = mc * dm; /* Leapfrog */ norm += pow(mc, -slope) * exp(-mc / mstar) * dmc - 0.5 * (pow(mc, -slope) / mstar + slope * pow(mc, -slope - 1.)) * exp(-mc / mstar) * (dmc * dmc); mc += dmc; } return norm; } /** * @brief Mean mass of Schechter function * * Expected cluster mass \int_{mmin}^{mmax} M N(M) dM */ double meanmass(const double mmin, const double mmax, const double slope, const double mstar, const double norm) { /* integration resolution*/ const int nintsteps = 10000; /* starting mass */ double mc = mmin; double meanm = 0.; /* Constant integration step in log space */ const double dm = pow(mmax / mmin, 1. / (double)nintsteps) - 1.; /* integrate to determine norm (total number) */ for (int n = 0; n < nintsteps; n++) { /* dm */ double dmc = mc * dm; /* Leapfrog */ meanm += pow(mc, -slope + 1.) * exp(-mc / mstar) * dmc - 0.5 * (pow(mc, -slope + 1.) / mstar + (slope - 1.) * pow(mc, -slope)) * exp(-mc / mstar) * (dmc * dmc); mc += dmc; } return meanm / norm; } /** * @brief Draw masses from Schechter function */ double schechtergen(const double mmin, const double mmax, const double slope, const double mstar, const double norm, const long long sid, const int cluster_index, const integertime_t ti_current) { /* integration resolution */ const int nintsteps = 10000; const double x = random_unit_interval_part_ID_and_index( sid, cluster_index, ti_current, random_number_mosaic_schechter); /* Method fails if x == 0 */ if (x == 0.) return mmin; /* starting mass */ double mc = mmin; double mcnorm = 0.; double dmc = 0.; double diff = 0.; /* Constant integration step in log space */ const double dm = pow(mmax / mmin, 1. / (double)nintsteps) - 1.; /* integrate until we reach the random number */ while (mcnorm / norm < x) { /* dm */ dmc = mc * dm; /* leapfrog */ diff = pow(mc, -slope) * exp(-mc / mstar) * dmc - 0.5 * (pow(mc, -slope) / mstar + slope * pow(mc, -slope - 1.)) * exp(-mc / mstar) * (dmc * dmc); mcnorm += diff; mc += dmc; } /* Interpolate between integration steps to correct mass */ /* Previous step */ const double x1 = (mcnorm - diff) / norm; /* Present step */ const double x2 = mcnorm / norm; const double frac = (x - x1) / (x2 - x1); /* interpolated mass */ return pow(mc - dmc, 1. - frac) * pow(mc, frac); } /** * @brief Draw masses from power-law */ double powerlawgen(const double mmin, const double mmax, const double slope, const long long sid, const int cluster_index, const integertime_t ti_current) { const double x = random_unit_interval_part_ID_and_index( sid, cluster_index, ti_current, random_number_mosaic_powerlaw); return pow(x * (pow(mmax, 1. - slope) - pow(mmin, 1. - slope)) + pow(mmin, 1. - slope), 1. / (1. - slope)); } /** * @brief Do star cluster evolution for the mosaics subgrid model * * @param sp The #spart to act upon. * @param e The #engine. * @param with_cosmology Are we running with cosmological time integration. */ void mosaics_clevo(struct spart* restrict sp, const struct engine* e, const struct cosmology* cosmo, const int with_cosmology) { const struct stars_props* props = e->stars_properties; const struct phys_const* phys_const = e->physical_constants; /* Proportional change in mass from stellar evolution */ const float stellar_evo_factor = sp->mass / sp->mass_prev_timestep; /* Total initial number of clusters over all model variations */ int ninit_sum = 0; /* Quick check to see if we need to bother going further */ for (int iMP = 0; iMP < MOSAICS_MULTIPHYS; iMP++) ninit_sum += sp->initial_num_clusters_evo[iMP]; if (ninit_sum == 0) { /* No clusters formed in this particle, just handle the field component */ for (int iMP = 0; iMP < MOSAICS_MULTIPHYS; iMP++) sp->field_mass[iMP] *= stellar_evo_factor; /* For stellar evo. mass loss at the next timestep */ sp->mass_prev_timestep = sp->mass; return; } /* Get the time since last step, then convert to physical units */ const integertime_t ti_step = get_integer_timestep(sp->time_bin); const integertime_t ti_begin = get_integer_time_begin(e->ti_current - 1, sp->time_bin); /* Calculate age of the star at current time */ double star_age; if (with_cosmology) { star_age = cosmology_get_delta_time_from_scale_factors( cosmo, (double)sp->birth_scale_factor, cosmo->a); } else { star_age = e->time - (double)sp->birth_time; } /* Get particle time-step */ double dt_star; if (with_cosmology) { dt_star = cosmology_get_delta_time(e->cosmology, ti_begin, ti_begin + ti_step); } else { dt_star = get_timestep(sp->time_bin, e->time_base); } /* Evaporation timescale constant */ double t0_evap; if (props->fixed_t0evap > 0) { /* Use a fixed value? */ t0_evap = props->fixed_t0evap; } else { /* Get t0 from local tidal field */ t0_evap = get_local_evaporation_t0(sp, props, cosmo->a3_inv); } /* Calculate the tidal heating parameters for this timestep */ double Mlheatsum[6] = {0.}; calc_heating(sp, star_age, dt_star, e->time, cosmo->a3_inv, Mlheatsum); /* Run multiple cluster formation models in parallel */ for (int iMP = 0; iMP < MOSAICS_MULTIPHYS; iMP++) { /* Initial number of clusters for this model, up to max array length */ const int ninit = min(sp->initial_num_clusters_evo[iMP], MOSAICS_MAX_CLUSTERS); if (ninit == 0) { /* No clusters were formed, just handle the field */ sp->field_mass[iMP] *= stellar_evo_factor; continue; } /* Any surviving clusters? */ if (sp->num_clusters[iMP] > 0) { /* Evolve the clusters */ for (int i = 0; i < ninit; i++) { if (sp->cl_mass[iMP][i] == 0.) continue; const double mass = (double)sp->cl_mass[iMP][i]; double dmevap = 0., dmsh = 0.; /* Mass loss by tidal shocks */ if (props->cluster_shocks) dmsh = mass_loss_shocks(sp, props, Mlheatsum, mass, phys_const->const_newton_G); /* Mass loss by evaporation */ if (props->cluster_evap) { const double initial_mass = (double)sp->cl_initial_mass[iMP][i]; /* Stellar particle metallicity */ const double logZ = log10(sp->chemistry_data.metal_mass_fraction_total + FLT_MIN); dmevap = mass_loss_evaporation( props, mass, initial_mass, t0_evap, dt_star, logZ, phys_const->const_newton_G, phys_const->const_solar_mass); } /* Total mass loss for this cluster */ double dmdis = dmevap + dmsh; /* Update cluster mass */ if (mass + dmdis < props->clMF_min) { /* Cluster went below minimum mass limit, consider fully disrupted */ /* fraction of mass loss */ dmevap *= -mass / dmdis; dmsh *= -mass / dmdis; /* new mass loss (can't lose more than it has) */ dmdis = -mass; /* We're done with this one */ sp->cl_mass[iMP][i] = 0.; sp->num_clusters[iMP]--; /* time of cluster disruption */ if (with_cosmology) { sp->cl_disruption_time[iMP][i] = cosmo->a; } else { sp->cl_disruption_time[iMP][i] = e->time; } } else { /* Cluster still survives, for now... */ sp->cl_mass[iMP][i] += (float)dmdis; } #ifdef SWIFT_DEBUG_CHECKS if (dmevap > 0 || dmsh > 0) { error("[mosaics] Positive mass loss?? dmevp = %g; dmsh = %g", dmevap, dmsh); } #endif /* Cluster mass lost to field */ sp->cl_dmshock[iMP][i] += (float)dmsh; sp->field_mass[iMP] += (float)fabs(dmdis); /* sp->cl_dmevap[iMP][i] += (float)dmevap; */ /* Stellar evolutionary mass loss for clusters */ sp->cl_mass[iMP][i] *= stellar_evo_factor; } /* Cluster loop */ } /* surviving clusters */ /* Apply stellar evolution to the 'field' populations */ sp->field_mass[iMP] *= stellar_evo_factor; for (int i = 0; i < ninit; i++) { sp->cl_dmshock[iMP][i] *= stellar_evo_factor; /* sp->cl_dmevap[iMP][i] *= stellar_evo_factor; */ } } /* End MOSAICS_MULTIPHYS loop */ /* For stellar evo. mass loss at the next timestep */ sp->mass_prev_timestep = sp->mass; } /** * @brief Do the cluster formation for the mosaics subgrid star cluster model * * @param sp The particle to act upon * @param stars_properties Properties of the stars model. * @param sf_props the star formation law properties to use * @param phys_const the physical constants in internal units. */ void mosaics_clform(struct spart* restrict sp, const struct engine* e, const struct cosmology* cosmo) { #if defined(STAR_FORMATION_NONE) error("Invalid SF model used!"); #endif const struct stars_props* props = e->stars_properties; const struct phys_const* phys_const = e->physical_constants; const double const_G = phys_const->const_newton_G; /* The resolved velocity dispersion (1/sqrt(3) converts to 1D assuming * isotropy) and add the subgrid thermal component. * Note: birth_velocity_dispersion is sigma^2 */ const double total_vel_disp2 = sp->sf_data.birth_velocity_dispersion * (1. / 3.) + sp->sound_speed_subgrid * sp->sound_speed_subgrid; /* -------- Get CFE -------- */ /* The local gas velocity dispersion (1D) */ float sigmaloc; if (props->subgrid_gas_vel_disp) { /* "Sub-particle turbulent velocity dispersion" (E-MOSAICS version) */ sigmaloc = sqrt(sp->birth_pressure / sp->birth_density); } else { /* The resolved velocity dispersion */ sigmaloc = sqrt(total_vel_disp2); } /* Sound speed of cold ISM. Fixed or subgrid value? */ float csloc; if (props->fixed_cs > 0) { csloc = props->fixed_cs; } else { csloc = sp->sound_speed_subgrid; } /* Stellar particle metallicity */ const double logZ = log10(sp->chemistry_data.metal_mass_fraction_total + FLT_MIN); /* Calculate the cluster formation efficiency based on local conditions * (Kruijssen 2012). */ const float CFE_local = f_cfelocal(sp->birth_density, sigmaloc, csloc, sp->cell_gas_fraction, logZ, props, phys_const); /* -------- Get Mcstar -------- */ /* The 'GMC' mass */ double M_collapse; if (props->Mgmc_Larson_relation) { /* Use inverted Larson-type M-sigma relation to get GMC mass */ M_collapse = props->GMC_m_sigma_const * pow(sigmaloc, props->GMC_m_sigma_index); } else { /* E-MOSAICS version: Toomre mass via tidal tensors (Pfeffer+18) */ /* Temporary array for eigvec/val calculation * Note the lower elements are never referenced by dsyevj3 */ double tensor[3][3]; tensor[0][0] = sp->tidal_tensor[2][0]; tensor[1][1] = sp->tidal_tensor[2][1]; tensor[2][2] = sp->tidal_tensor[2][2]; tensor[0][1] = sp->tidal_tensor[2][3]; tensor[0][2] = sp->tidal_tensor[2][4]; tensor[1][2] = sp->tidal_tensor[2][5]; /* Get the eigenvalues */ double tideval[3]; dsyevj3(tensor, tideval); /* sort eigenvalues, tideval[2] (=lambda_1) is largest eigenvalue */ sort3(tideval); /* Circular and epicyclic frequencies */ double Omega2; if (props->Omega_is_lambda2 && tideval[2] < 0) { /* Correct version of Omega for a spherical distribution, in principle. * When tideval[2] > 0 we can't guarantee that tideval[1] will be < 0 */ Omega2 = fabs(-tideval[1]); } else { /* For axisymmetric system, derived in Pfeffer+18 (i.e. E-MOSAICS version) */ Omega2 = fabs(-tideval[0] - tideval[1] - tideval[2]) * (1. / 3.); } double kappa2 = fabs(3. * Omega2 - tideval[2]); /* Factor out cosmological factors from birth values */ Omega2 *= cosmo->a3_inv; kappa2 *= cosmo->a3_inv; sp->Omega_birth = sqrt(Omega2); sp->kappa_birth = sqrt(kappa2); /* Pressure from turbulent and thermal contributions */ const double total_pressure = sp->birth_density * total_vel_disp2; /* Gas surface density, assuming hydrostatic equilibrium */ const double surfDens = sqrt(2. * sp->cell_gas_fraction * total_pressure / (M_PI * const_G)); /* The Toomre mass */ const double MTconst = 4. * M_PI * M_PI * M_PI * M_PI * M_PI * const_G * const_G; sp->Toomre_mass = MTconst * surfDens * surfDens * surfDens / (kappa2 * kappa2); /* Total collapse time of Toomre volume */ const double tcollapse = sqrt(2. * M_PI / kappa2); /* Feedback timescale (Kruijssen 2012) */ double tfb = feedback_timescale(sp->birth_density, sigmaloc, csloc, props, phys_const); /* Toomre collapse fraction (Reina-Campos & Kruijssen 2017) */ sp->frac_collapse = fmin(1., tfb / tcollapse); sp->frac_collapse *= sp->frac_collapse; /* f^2 */ sp->frac_collapse *= sp->frac_collapse; /* f^4 */ M_collapse = sp->Toomre_mass * sp->frac_collapse; } /* ----- Set up the star cluster population for this particle ----- */ /* Run multiple cluster formation models in parallel */ for (int iMP = 0; iMP < MOSAICS_MULTIPHYS; iMP++) { /* Initialise */ sp->initial_num_clusters[iMP] = 0; sp->initial_num_clusters_evo[iMP] = 0; sp->num_clusters[iMP] = 0; sp->initial_cluster_mass_total[iMP] = 0.; sp->initial_cluster_mass_evo[iMP] = 0.; for (int i = 0; i < MOSAICS_MAX_CLUSTERS; i++) { sp->cl_mass[iMP][i] = 0.; sp->cl_initial_mass[iMP][i] = 0.; sp->cl_dmshock[iMP][i] = 0.; sp->cl_disruption_time[iMP][i] = -1.; /* sp->cl_id[iMP][i] = -1; */ /* sp->cl_dmevap[iMP][i] = 0.f; */ } /* The cluster formation efficiency for this formation model */ if (props->fixed_CFE[iMP] > 0) sp->CFE[iMP] = props->fixed_CFE[iMP]; else sp->CFE[iMP] = CFE_local; /* Exponential truncation of mass function for this formation model */ if (props->power_law_clMF[iMP]) { /* No truncation for pure power law */ sp->Mcstar[iMP] = 0.; } else { /* Fixed value or locally varying? */ if (props->fixed_Mcstar[iMP] > 0) sp->Mcstar[iMP] = props->fixed_Mcstar[iMP]; else sp->Mcstar[iMP] = props->GMC_SFE * sp->CFE[iMP] * M_collapse; } /* Power-law index of the mass function, M^(-alpha) */ const float clMF_slope = props->clMF_slope[iMP]; /* Cluster and field masses if we were conserving mass within particles */ const float part_clust_mass = sp->mass_init * sp->CFE[iMP]; sp->field_mass[iMP] = sp->mass_init - part_clust_mass; /* Determine mean cluster mass by integrating mass function, * along with the integral normalisation */ double mean_mass = 0.; double norm = 0.; if (props->power_law_clMF[iMP]) { /* Power law. Treat some special cases first */ if (clMF_slope == 1.) { norm = log(props->clMF_max) - log(props->clMF_min); mean_mass = props->clMF_max - props->clMF_min; } else if (clMF_slope == 2.) { norm = 1. / props->clMF_min - 1. / props->clMF_max; mean_mass = log(props->clMF_max) - log(props->clMF_min); } else { norm = (pow(props->clMF_max, 1. - clMF_slope) - pow(props->clMF_min, 1. - clMF_slope)) / (1. - clMF_slope); mean_mass = (pow(props->clMF_min, 2. - clMF_slope) - pow(props->clMF_max, 2. - clMF_slope)) / (clMF_slope - 2.); } mean_mass /= norm; } else { /* Schechter cluster mass function */ /* If Mcstar < clMFmin we have undefined behaviour and just form all * clusters at M=clMFmin */ if (sp->Mcstar[iMP] > props->clMF_min) { /* Get the normalisation of the integral */ norm = schechternorm(props->clMF_min, props->clMF_max, clMF_slope, sp->Mcstar[iMP]); } if (norm > 0) { mean_mass = meanmass(props->clMF_min, props->clMF_max, clMF_slope, sp->Mcstar[iMP], norm); } } /* Number of clusters to try and form. * Draw from Poisson distribution, so we can be resolution independent */ if (mean_mass > 0) { const double lambda = round(part_clust_mass / mean_mass); sp->initial_num_clusters[iMP] = random_poisson( sp->id, lambda, e->ti_current, random_number_mosaic_poisson); } else { sp->initial_num_clusters[iMP] = 0; } /* Lowest mass cluster we've formed so far */ float min_formed_mass = FLT_MAX; int min_mass_idx = -1; /* Counter for our cluster arrays */ int iArr = 0; /* draw cluster masses from mass function */ for (int i = 0; i < sp->initial_num_clusters[iMP]; i++) { /* Try and form a cluster */ float mTry; if (props->power_law_clMF[iMP]) { /* Power law mass function */ mTry = powerlawgen(props->clMF_min, props->clMF_max, clMF_slope, sp->id, i, e->ti_current); } else { /* Schechter cluster mass function */ mTry = schechtergen(props->clMF_min, props->clMF_max, clMF_slope, sp->Mcstar[iMP], norm, sp->id, i, e->ti_current); } sp->initial_cluster_mass_total[iMP] += mTry; #ifdef SWIFT_DEBUG_CHECKS if (!isfinite(mTry)) { error("[mosaics] mTry is not finite!! mTry=%g, Mcstar=%g, norm=%g", mTry, sp->Mcstar[iMP], norm); } #endif /* Lower than our limit to evolve clusters */ if (mTry < props->clMF_min_evolve) { /* Add destroyed clusters to field */ sp->field_mass[iMP] += mTry; continue; } sp->initial_cluster_mass_evo[iMP] += mTry; /* Check if we've run out of space... */ if (iArr < MOSAICS_MAX_CLUSTERS) { /* Still space left */ /* sp->cl_id[iMP][iArr] = iArr; */ sp->cl_mass[iMP][iArr] = mTry; sp->cl_initial_mass[iMP][iArr] = mTry; /* Lowest mass cluster we've formed so far */ if (mTry < min_formed_mass) { min_formed_mass = mTry; min_mass_idx = iArr; } } else { /* Ok we've run out of space. Replace the least massive * cluster in the array so we at least track the most massive ones */ if (mTry > min_formed_mass) { /* Add the previous one to field */ sp->field_mass[iMP] += sp->cl_mass[iMP][min_mass_idx]; sp->cl_mass[iMP][min_mass_idx] = mTry; sp->cl_initial_mass[iMP][min_mass_idx] = mTry; /* reset the counters */ min_formed_mass = FLT_MAX; min_mass_idx = -1; for (int j = 0; j < MOSAICS_MAX_CLUSTERS; j++) { if (sp->cl_mass[iMP][j] < min_formed_mass) { min_formed_mass = sp->cl_mass[iMP][j]; min_mass_idx = j; } } } else { /* We're out of space... just add to field */ sp->field_mass[iMP] += mTry; } } iArr++; } /* Current number of clusters */ sp->num_clusters[iMP] = min(iArr, MOSAICS_MAX_CLUSTERS); /* Number of clusters above the evolution mass limit */ sp->initial_num_clusters_evo[iMP] = iArr; /* Warn when we hit the cluster array end */ if (sp->initial_num_clusters_evo[iMP] > MOSAICS_MAX_CLUSTERS) { message( "sp->id=%lld: Trying to form more clusters (%d) than allowed (%d)\n" "MULTIPHYS=%d; mass_init=%g; clust_mass=%g; mean_mass=%g", sp->id, sp->initial_num_clusters_evo[iMP], MOSAICS_MAX_CLUSTERS, iMP, sp->mass_init, part_clust_mass, mean_mass); } } /* End MOSAICS_MULTIPHYS loop */ }